{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Singularities due to Wood's Anomaly\n",
    "Wood's anomly is an age-old phenomenon observed in the early 20th century. It is basically a phenomenon when $k_z = 0$ in a layer. \n",
    "In RCWA, we define a $k_{zi}$ for every Fourier component:\n",
    "$$\n",
    "k_{zi} = k_0^2n^2 - k_{xi}^2 -k_{yi}^2\n",
    "$$\n",
    "Assume for a second we are in the 1D case and $k_y = 0$ and that $n=1$.Then we can make $k_{zi} = 0$ simply by asking $k_0 = k_{xi}$, or:\n",
    "$$\n",
    "\\frac{2\\pi}{\\lambda} = k_x \\pm \\frac{2\\pi m}{a}\n",
    "$$\n",
    "\n",
    "Further assuming normal incidence, we just have:\n",
    "$$\n",
    "\\frac{2\\pi}{\\lambda} =  \\frac{2\\pi m}{a}\n",
    "$$\n",
    "\n",
    "So we see we can get issues precisely when $\\lambda$ is a rational fraction of the lattice constant. However, for a typical grating, it's not really obvious how we can get to this singularity because the $k_z$ values we are interested are extracted from an eigensolver. One case we can force $k_z=0$ without a doubt it for a uniform slab. However, as the below script will show, the 1D TE and TM Case with the Gaylord formulation seems to be safe, specifically because unlike the scattering matrix formalism, they do not try to invert $Kz$ when they go extract the eigenmodes in $H$\n",
    "\n",
    "This is the link to the original paper\n",
    "https://www.osapublishing.org/view_article.cfm?gotourl=https%3A%2F%2Fwww%2Eosapublishing%2Eorg%2FDirectPDFAccess%2FE451A175-DDA2-A95A-D5C67C8D01CB3352_33172%2Fjosaa-12-5-1068%2Epdf%3Fda%3D1%26id%3D33172%26seq%3D0%26mobile%3Dno&org=Stanford%20University%20Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd809X6wPHPSZruFrpYLaMs2bMi\niCgqoqIX90DcA7dXvQ5cV67+9F73Xrg3Km5AFJUhImhBQJBRWlaZXXSvJOf3xzdJ0zZt0jZt2vK8\nXy9eLek33+9pKE9Pnu9znqO01gghhGg/TIEegBBCCP+SwC6EEO2MBHYhhGhnJLALIUQ7I4FdCCHa\nGQnsQgjRzkhgF0KIdkYCuxBCtDMS2IUQop0JCsRF4+Pjda9evQJxaSGEaLNWr16drbVO8HZcQAJ7\nr169SE1NDcSlhRCizVJK7fTlOEnFCCFEOyOBXQgh2hkJ7EII0c4EJMcuhDj8VFZWkpmZSVlZWaCH\n0uqFhoaSlJSExWJp1PMlsAshWkRmZiZRUVH06tULpVSgh9Nqaa3JyckhMzOT5OTkRp1DUjFCiBZR\nVlZGXFycBHUvlFLExcU16Z2NBHYhRIuRoO6bpr5OEtiFqMPqvzby+QcvU2mzB3ooQjSIBHYh6rDn\nu6c4Z9s9rNni05oQ0QY8//zzDBw4kGnTpjFp0iRGjBjBJ598wrPPPktJSUmgh+c3cvNUiDoklBsB\nvezAVhjUuJtYonV5+eWX+e677zhw4AB33303a9euBYzV8BdffDHh4eEBHqF/SGAXog6Jtt0AqJw0\n4OTADkY02XXXXUdGRgZTpkxh69atREZGMmLECK644gr27t3L8ccfT3x8PIsXLw70UJtMArsQHljL\nS0jUB0FB8KH0QA+n3ek1c36znHfH/06r82uvvvoqCxcu5JdffmHDhg08+eSTzJs3D4BnnnmGxYsX\nEx8f3yzjammSYxfCg8K9WzErDUBkYUaARyNEw8iMXQgPyvZtAiBTxxNbuiOwg2mH6ptZi6aTGbsQ\nHlgPpgHwk20kCRV7wG4L8IhEc4qKiqKwsDDQw/AbCexCeGDKTWOPjmODTsZCJRySksf2bMaMGZx6\n6qkcf/zxgR6KX0gqRggPQvPT2WjvRrq9m/FAdhrE9g7soEST7dixA4CJEycyceJE1+M333wzN998\nc2AG1Qxkxi5ETVoTVbSddN2NdG0Edp21JcCDEsJ3EtiFqKlwH8G2EtJ1N/KJJEtHU3lAArtoOySw\nC1FT9lYA12w9Q3fDnrU1kCMSokH8kmNXSu0ACgEbYNVap/jjvEIERLZREePMr6fbuzIqb10gRyRE\ng/jz5unxWutsP55PiMDITqOEMA7SkdiIYNLLumEpWwwluRAeG+jRCeGVpGKEqCl7KztVIqDokxDh\nSsk4Z/JCtHb+Cuwa+EEptVopNcNP5xQiMLLTSLN1BaBPQqRbYJc8uzC888473HTTTV6P2bt3r+vv\nV199NX///XdzDw3wXypmvNZ6r1KqE7BIKbVZa73M/QBHwJ8B0KNHDz9dVgg/Ky+Cgky22I4m2Gyi\ne2w4e3QCVhVMkAR20QDvvPMOQ4YMoVs3Y2LwxhtvtNi1/TJj11rvdXw8CHwJjPFwzGytdYrWOiUh\nIcEflxXC/3K2AUZFTEyEhbiIYOyYyApOklRMO3HmmWcyevRoBg8ezOzZswGIjIzkvvvuY/jw4Ywd\nO5YDBw4A8O2333LUUUcxcuRIJk2a5HrcqbCwkOTkZCorKwEoKCigV69efPbZZ6SmpjJ9+nRGjBhB\naWkpEydOJDU1FYCFCxcyatQohg8fzoknnuj377HJM3alVARg0loXOj6fDDzU5JEJEQjOihjdjZjw\nYGIiggHYY06iq8zY/ee7mbD/L/+es8tQOPV/Xg976623iI2NpbS0lCOPPJJzzjmH4uJixo4dyyOP\nPMJdd93F66+/zv33388xxxzDypUrUUrxxhtv8Pjjj/PUU0+5zhUVFcXEiROZP38+Z555JnPmzOGc\nc87hvPPO46WXXuLJJ58kJaV6kWBWVhbXXHMNy5YtIzk5mdzcXP++Dvhnxt4ZWK6UWgf8DszXWi/0\nw3mFaHnZW9HKxE7dmdiIYGIdgX073SBvB1grAjs+0WTPP/+8a2a+e/du0tLSCA4O5vTTTwdg9OjR\nrtYDmZmZnHzyyQwdOpQnnniCjRs31jrf1Vdfzdtvvw3A22+/zRVXXFHv9VeuXMmxxx5LcrKxK1ds\nrP8rrZo8Y9daZwDD/TAWIQIvJ43i8CQqSi3VAvtWW1fQNsjbDglHBHiQ7YAPM+vmsGTJEn788Ud+\n++03wsPDmThxImVlZVgsFpRSAJjNZqxWK2D0kLn99tuZOnUqS5YsYdasWbXOOX78eHbs2MHSpUux\n2WwMGTKk3jForV3Xai5S7iiEu+w08sJ6AhiBPdwI7OvLOzu+LumYtiw/P5+YmBjCw8PZvHkzK1eu\n9Hp8YmIiAO+++26dx1166aVMmzat2my9rlbA48aNY+nSpWzfvh2g1aZihGgf7DbI2caBYKNqKyY8\nmOgwCyYFf5V1Mo6RwN6mnXLKKVitVoYNG8YDDzzA2LFj6z1+1qxZnHfeeUyYMKHebfOmT59OXl4e\n06ZNcz12+eWXc91117lunjolJCQwe/Zszj77bIYPH84FF1zQ9G+sBqW19vtJvUlJSdHOu8NCtBp5\nO+C54XyZdDe3bRvOQ2cM5tJxvRj18CJyiytIj78Dc5/j4KxXAz3SNmnTpk0MHDgw0MNoFnPnzuXr\nr7/m/fff99s5Pb1eSqnVvrRskX7sQjg5KmIyMOqOYxxpmNiIYHKLKyjr0IcImbGLGm6++Wa+++47\nFixYEOihuEhgF8LJEbQ3VXYBbK4bp848e2FkLyJ2fA1aQzPf/BJtxwsvvBDoIdQiOXYhnLLTICyW\nXWVhQNWMPSbCAkBOaE8oL4CigwEbYlsXiNRvW9TU10kCuxBO2WkQ35/cYmMVYVykMxUTAsA+i6MV\nhqRjGiU0NJScnBwJ7l5orcnJySE0NLTR55BUjBBO2VvR/SaTl24sQuoYbszUYx0z9l2mRNdxJE8I\nyBDbsqSkJDIzM8nKygr0UFq90NBQkpKSGv18CexCAJTmQfFByjr0wWbXRIYEERJkBqpSMrsqO4Il\nQnrGNJLFYnGtthTNS1IxQgBkG82/8iMdy7wdN07dP88rrYT4fpKKEa2eBHYhwBWss0OMVacxHgJ7\nbnGFI7DLjF20bhLYhQAjsJss7DcZrQNiHfl1qBnY+0P+LqgoCcgwhfCFBHYhwOjDHteH3FI7UH3G\n7syx5zln7AC56S0+RCF8JYFdCDBm7PH9yC0xKmLiPKViShwzdufxQrRSEtiFsFVCbgbE9TNm5VSf\nsYcHmwkJMlFWaackqiegJM8uWjUJ7ELk7QC7FeL7k+MI7M42AgBKqapZe7kJYnrKjF20ahLYhXAG\n6fj+rhm7e7kjuOfZKyFOSh5F6yaBXQhXYO/ryrHXDOzO9gI5xeVGnj17G9jtLTpMIXwlgV2I7DSI\n7AKhHTzm2MFtxl7iqIyxlkJBZosPVQhfSGAXIjvNVcboKccO7rXslW6VMXIDVbROEtjF4U1rR6lj\nfyptdgrLrJgUdAizVDusei27BHbRuklgF4e34mwoO2TcOHXk12PCgzGZqm+kEevKsVdARDyEdpQb\nqKLV8ltgV0qZlVJ/KqXm+eucQjQ79xundeTXoSo1k1dcYeyeFN9fArtotfw5Y/8nsMmP5xOi+bmV\nOubWUeoIVbsoOatmpBmYaM38EtiVUknAacAb/jifEC0mOw2CwiA6yahRp/aNU6jRCAyMwF60H8ry\nW2yoQvjKXzP2Z4G7gDoLe5VSM5RSqUqpVNlBRbQa2Vshvi+YTK7ZuMdUTIRbKgbcbqBua5FhCtEQ\nTQ7sSqnTgYNa69X1Hae1nq21TtFapyQkJDT1skL4R06aK0jnFjlTMZZah7nXsdvtuiqw50g6RrQ+\n/pixjwemKqV2AHOAE5RSH/jhvEI0r8oyyNvpCtJ5rlWnIbUOtZhNRIUGYddQUFYJMb3AFCQ3UEWr\n1OTArrW+R2udpLXuBVwI/Ky1vrjJIxOiueWmA9q1OKnq5mntGTtUtfLNKa4AswVie0tgF62S1LGL\nw1fWFuNjnBHY3evYPYnxlGfPksAuWh+/Bnat9RKt9en+PKcQzSZ7K6Cq2gkUOTfZqJ2KgapqmWqV\nMbnpRj93IVoRmbGLw1fWFujYAyxhgNuMvY5UTEytkscjjD7uudubf6xCNIAEdnH4yt4KCUcAoLWu\nd4ESVOXYXYuUEpwlj1uad5xCNJAEdnF4stuMDawdFTElFTbKrXZCgkyEWcwen+Ixxw5yA1W0OhLY\nxeHp0C6wlrlm7M7ZelxEMEopj0+pyrE7cuohURCdKDdQRasjgV0cnlw9YozAnlfPqlOnqrYC5VUP\nxveXVIxodSSwi8OTs9TRkSf3ll8Ht5unJW5VMAlHGP1mtG6ecQrRCBLYxeEpewtEdIKwGMC3wF6r\nXwwYM/aKIijY03xjFaKBJLCLw1NWVUUMVAX2uhYngYc6dqg6R5akY0TrIYFdHH60NmbszqoW3PvE\n1B3Yo8OCMJsUReVWyq0240FHjl4qY0RrIoFdHH6KDhp91D3N2OsJ7Eop14z+kDPPHhFvpHNkxi5a\nEQns4vDjrGJxm7G7lzvWx9kgzJWOUcqYtcuMXbQiEtjF4cdVEVM1Y3funlRfjt3969Xz7P1lxi5a\nFQns4vCTvRWCoyCqq+uhXB9y7ABxkR4Ce3x/KMmGklz/j1WIRpDALg4/WVuMzoxuK0x9KXeE6jsp\nucRLZYxoXSSwi8NPdvVSR5tdc8gRqDuGe+7s6OQM/M4Wv4A0AxOtjgR2cXgpy4fCfdVunBaUVmLX\nEB0ahMVc/38J1yIl9xl7hx4QFCY9Y0SrIYFdHF6yHZtPu5c6+phfdz+mWo7dZIL4vjJjF62GBHZx\neHHmweNr17D7Etg95tid55MZu2glJLCLw0v2FjAHQ0wv10MNCewec+xgvAPI3wUVxX4bqhCNJYFd\nHF6ytkJsHzAHuR7K86FPjJPHHDu4bbqR5p9xCtEEEtjF4SV7S1UVi0OOc8Ye2YBUTHEl2r1Vb4L0\njBGtR5MDu1IqVCn1u1JqnVJqo1LqP/4YmBB+V1kGeTuq5dehasYe68OMPSzYTJjFTIXNTlG5teoL\nsX1AmaWWXbQK/pixlwMnaK2HAyOAU5RSY/1wXiH8KzcdtL1aRQxUVcXU1wDMXVVfdrcNN4KCITZZ\nKmNEq9DkwK4NRY6/Whx/ZDsZ0fpk1W7+BQ2bsYNbyaNUxohWyi85dqWUWSm1FjgILNJar/LHeYXw\nq+ytgDLaCbjJbUCOHapm9tV2UgIjd5+bDrZKD88SouX4JbBrrW1a6xFAEjBGKTWk5jFKqRlKqVSl\nVGpWVpY/LitEw2RtgY49wBJW7WHXAiVfZ+yOtgM5NQN7/BFgt0Lu9qaPVYgm8GtVjNb6ELAEOMXD\n12ZrrVO01ikJCQn+vKwQvqnRI8bJ1bLX5xx7iON5HmbsIHl2EXD+qIpJUEp1dHweBkwCNjf1vEL4\nld1m1JjXyK+XW20UlVsJMimiQ4PqeHJ1rs026qpll8oYEWC+/STXryvwrlLKjPGL4lOt9Tw/nFcI\n/zm0E2zltWbs7rN15dbGtz7OmX1uzdWnIVEQnSi17CLgmhzYtdbrgZF+GIsQzcdZrVKjhj23gRUx\n7sfWmrGDMWuXGbsIMFl5Kg4Pzrx3jVWneQ3o7OgUW1dVDBjvCLLTwG5v3DiF8AMJ7OLwkLUVIjpB\nWEy1h3Ma0ADMqc46djBm7JXFULCn8WMVookksIvDQ/aWOipinKtO6985yV2Mp57sTq6eMZKOEYEj\ngV20f1obM/YaFTHQuBx7xzALSkF+aSVWW42Ui2v/U7mBKgJHArto/4oOQHm+5xl7I3LsQWYTHcIs\naG0E92oi4o10j8zYRQBJYBftXx09YqAqx+7r4iQnV2VMzXSMUtIzRgScBHbR/jnryuvJsTdkxg7e\n8uz9ZcYuAkoCu2j/srdCcBREda31pYZsi+euzp2UwHhnUJIDxTkNH6sQfiCBXbR/WVuMjo4eVpY2\nOrA7UjG1GoFB1Q1UWYEqAkQCu2j/6mj+pbV2zbh92e/UXZ2te0GagYmAk8Au2reyfCjc5/HGaVG5\nlUqbJjzYTKjF3KDTuhqBFXvovd6hBwSFyQ1UETAS2EX7lp1mfKynXW9D0zDGcxytez3l2E0miO8r\nM3YRMBLYRfvmKnWsHdhzisuBxgb2OjbbcJKSRxFAEthF+5a9BczBENOr1pcam193f47HHDsY7xDy\nd0FFcYPPLURTSWAX7VvWVojtA+baHaqd+fG4RszY4xypGI917FCV03emgoRoQRLYRfuWvaVWq16n\nXEcqpqGrTo3nOG+e1jNjByl5FAEhgV20X5VlkLfDY34dqmbsMeG+d3Z0igwJwmJWlFbaKK2w1T4g\ntg8os2y6IQJCArtov/K2g7Z7LHWEqhl7fGRIg0+tlKpKx3iqjAly5PVzJBUjWp4EdtF+5aQbH+N6\ne/6yY8/SuEYEduN5jtWnReV1HNAHcjIadW4hmkICu2i/ch1BNdZzYM8udgb2hufYjecZvxByam5q\n7RTbxxiD1o06vxCNJYFdtF+56RAWW2s7PCfnTDs+onEz9njHTdfsumbssb2NbfKKDjTq/EI0lgR2\n0X7lpBvpkLq+XNTUGXs9jcCgKgXkTAkJ0UKaHNiVUt2VUouVUpuUUhuVUv/0x8CEaLLc7UY6xIOS\nCiullTZCLSbCgxvWJ8apKhVT14zdce1cybOLllV71UbDWYF/aa3XKKWigNVKqUVa67/9cG4hGqey\nFAoy68yvu2brESEoD+18feFsRVBnjr1DdzAFGSkhIVpQk2fsWut9Wus1js8LgU1AYlPPK0ST5G43\nPtaRinHmxRubhgGIdzw3u65UjDnIUfIogV20LL/m2JVSvYCRwCoPX5uhlEpVSqVmZWX587JC1Oal\nIqZqxt74wO6sY68zFQOOypjtjb6GEI3ht8CulIoEPgdu1VoX1Py61nq21jpFa52SkJDgr8sK4Zkz\n/VFXYC92ztgbVxFjPNdLKsZ5fSl5FC3ML4FdKWXBCOofaq2/8Mc5hWiSnHQIj4Owjh6/nN3Eihhw\nm7EXl6PrCtxxfYySx8L9jb6OEA3lj6oYBbwJbNJaP930IQnhB7kZdVbEQNUsu7E17ABhwWYigs1U\n2jQFZVbPBznfMUhljGhB/pixjwcuAU5QSq11/Jnih/MK0Xi5GXWmYcA9FdP4GbvxfG8lj87ALjdQ\nRctpcrmj1no50Lh6MSGaQ0UJFOzxcXFS42fsxvOD2ZVbQk5xBb093Trq0B1MFqmMES1KVp6K9idv\nh/Gxnhm7q9yxCVUxxvO9zNidJY+SihEtSAK7aH+8VMRAVRuAxrTsdeeqZfelMkaIFiKBXbQ/rna9\nnlMxdrt27XzUmI2s3flU8hgnXR5Fy5LALtqf3HQIj4fQDh6/nF9aic2uiQ4NIjioaf8F3Ese6xTb\nGypLoHBfk64lhK8ksIv2J3d7/TdOm7BzUk0+z9hB0jGixUhgF+1PTrqXG6dNX5zk5PzlUGdPdqga\ni1TGiBYigV20LxUlULjXp8VJcU1YnOTktSc7GCWP5mCpZRctRgK7aF/ynF0dm39xEvjYCMxklpJH\n0aIksIv2Jcd7qWO2nxYnAcSEW1AK8koqsdrsdR8Y21s2thYtRgK7aF9cNez1pWKcN0+bPmMPMpuI\nCTfOk1tSXy27o+TRXk/wF8JPJLCL9iU3AyISIDS6zkP8mWM3zuNLZUxvsJZCkXR5FM1PArtoX3Lq\nb/4F/s2xu5/Ha192kMoY0SIksIv2JTe93jQMuLXs9Vtg92WRkrOWXQK7aH4S2EX7UVFsrO6spyIG\n3BuA+ScVEx/hQ7+YDkmOkke5gSqanwR20X449xatJxVTYbVTUGbFbFJ0CLP45bJee7JDVcmjpGJE\nC5DALtoPHypi3Jt/mUz+2UbApxy7c1wyYxctQAK7aD+cQbMF+rC786kRGDi6PG6XkkfR7CSwi/Yj\nJ917qaOf+rC786knO0BsslHyKF0eRTOTwC7aDy8bWENVHtxfpY7GuXycsUtljGghEthF+5GbUW+7\nXvD/4iRoQI5d2veKFiKBXbQPzlJHL4uTsv28OAkgKiSIYLOJkgobJRXWug+MTgJziFTGiGbnl8Cu\nlHpLKXVQKbXBH+cTosF8uHEK/l+cBKCU8m3WbjJJl0fRIvw1Y38HOMVP5xKi4ZzB0msqxr+Lk5x8\n6ssOVfufCtGM/BLYtdbLgFx/nEuIRvGhXS9UBV5/pmLAx77sYIxPujyKZiY5dtE+5KZDRCcIiar3\nsKpUTDPN2L2WPPYGa5mxy5MQzaTFArtSaoZSKlUplZqVldVSlxWHCy8bWANorasWKPl5xu7a+9SX\nRUog6RjRrFossGutZ2utU7TWKQkJCS11WXG4yPHe1bG4wka51U6YxUx4cJBfL+9TT3aoGqNUxohm\nJKkY0faVFxkbWMQm13tYcyxOcvKpERhAdKJR8iiLlEQz8le548fAb8ARSqlMpdRV/jivED5xbWBd\n/4zdn3ud1uRzVYzJZPwCcnaiFKIZ+OX9qNZ6mj/OI0Sj5Hjv6ghue536sQGYU7yjKsZrvxgwximp\nGNGMJBUj2j5Xu14vqZhmKnV0P6fXVAwY48yTLo+i+UhgF21fbgZEdvah1NGZY/d/KibW8S4gt7gC\nu13Xf3BcHyl5FM1KArto+3K8d3UEtxx7M6RiQi1mokKCsNo1BWWV9R8slTGimUlgF21fbrrXFafQ\nPL3Y3cX53JfdMVapjBHNRAK7aNvKi6DogNcNrKF5yx2N8zag5DEoVBYpiWYjgV20ba6ujt5TMc3R\ni92da5GSLyWPMclGCkmIZiCBXbRtznSGlxp2qNrhyJ8te935PGMHR5dHScWI5iGBXdTPZoWNXwZs\nQc2jCzbx4Ncb0LqOShPnDciY+ksdbXZNrmMmHdMMN0+hAXufQtUipTpKHtdnHuKi11eyeX+BP4fo\nG61hwxdyc7cNk8AuPNMaNs+HV46Gzy6HD881dilqQat35jJ7WQbv/raT7dl1XDt3O0R2gZDIes91\nqKQCu4aO4RYs5ub5sa9KxfhSy94HbOVQsMfjl9/4ZTsr0nN4ZP4mfw7RN2veg7lXwEtjYP4dUCRN\n+9oaCeyitt2/w9unwpyLQNvhhAeM2dvCmS06jNeWVuWgV22vo91/brqPaZjmK3V0qkrF+DBjj6t7\nY2utNau25wDwS1o2G/fm+22MXmWnGf/OvSbAqMsg9S14fgQsecy4US3aBAnsokp2GnxyMbx5knFT\n8vRn4IaVcOwdcMytxkzu769bZCjbDhaxaNMB199/9xTYtTbG7GXFKeDWrrd5bpwa5/axwyNU3ezN\nTqv1pV25JRwoqJr1z17WQjdZrRXw+dUQFAJnz4bTn4Ybf4c+J8CSR+H5kfDHm2DzUqcvAk4CuzDe\nas+7DV46CtIXw/H3wS1/QsqVYHa0E5p4L3QbCd/cAvme0wf+9MYvGWgNo3vGALAqI6d2nj1/N5Rk\nQ9cRXs/XHHud1uRzT3aA6G4QHg97/6z1pVUZxi+x4UkdMJsU89bvIzOvxK9j9WjxI7BvLUx9wRgf\nQHxfuOB9uGqR8S5j/u3w8ljY9G3zj0c0mgT2w52tEt47w5iNH3kV3LIWjrsLgiOqHxcUDOe8CbYK\n+Oq6Zu1zcrCwjC/W7EEpeOycYXQIs7A3v4zMvNLqB+5ZbXxMHO31nM2116k7n3uyAyhljNv5PbhZ\n6UjDTB2RyD+GdcVm17y5vJlvXm9fBr8+Z6RfBv6j9te7j4ErvoMLPwZlNt7Zbfi8ecckGk0C++Hu\njzfg4EY47x2Y8gRE1rMJSlwfOPUxIwiseL7ZhvTOrzuosNmZPKgzfTtFcmSvWMBDnj0z1eht3nmI\n13M2ZwMwp47hwZgU5JdWUmnz4RdfUgpkbYGy6pUvzrTTUcmxzDjWSNnM+X03ed7q4xurJBe+uBbi\n+sIp/637OKVgwBS4fgV0HQ7f3y9591ZKAvvhrOggLH4U+k6CAaf79pyRl8DAqfDzwx7TCE0eUrmV\n91fuBHAFtbG9HYE9I6f6wXtWGwEmyHuwbs5e7E5mk3I1A/MpCCeOBjTsXeN6aM+hUjLzSokKDWJg\n12gGdYtmQr94SittfOB4XfxKa/j2FijOgnPeqP1OzRNzEEx50mhi9suT/h+TaDIJ7P5Qmgc/3A8f\nXQhbFraddqw/zoLKUjjlMWM25gul4B/PGd0UP7/a7yWQc37fRWGZlSN7xbjy60clxwE1Zuw2K+xd\n61MaBpq3F7u7uIb0ZU8cZXx0S8c4f3mN6RWL2WT8m1x3nPEL7t3fdlBWafPfYMFIwW36Fk58ALp5\nv1fh0n0MjJgOK16E7G3+HVNzKTwAP/4HPjgXtv0U6NE0KwnsTWGrhJWvGtUCK140Zl4fX2DUfq/9\nyKgyaK12/w5rP4SjbzJukDVEeCyc9ZrfSyArbXZXLvnaY6tKGAd2jSIyJIhduSXsy3fk2Q/+DdZS\nI53hg6pUTPPN2I3zN6CWPSzGSH9kVgV2ZxpmTHKs67Gj+8QxJDGa7KIKPl+T6b/BOksbk4+FcTc3\n/PmTZoElDL67y5j5t1bZ24yb/s8OgV+fhX3r4IOz4YNz4GAA1gm0AAnsjeFcvPPyWFh4t5EOuO4X\nuG0jnP06mMzw1fVG/e+KF6G8MNAjrs5ugwV3QFQ3mHBH486RPMHvJZDfrtvLvvwy+naK5IQBnVyP\nB5lNpPQyZu+ussc9qcbHBs7YmzPHbpy/AbXsAIkpxvfiCIzOdyVH9Y5zHaKUcv2ie31ZBjZv/d59\n4V7aeNZrRv+ahorsBMffC+k/wZYFTR+Tv2Wmwpzp8GIKrP/ESCPelAq3bYDJ/we7/zAmYfNua3eL\nsPy7VfvhYO9aI+2y4xeI7w8XfQb9TqpKZQw7H4aeZ7zV+/VZ+OE+WPq4UXEy7kaIiPf7kD5ctZPX\nlmbwr8n9OWNEovcnrHnXmLUeBvT/AAAgAElEQVSc+5bXFZsAWYXl7M8vo8xqo7zSTlmljXKrnfIO\nl3Jchx+I/PImSmKGE9O1V6O/B621q157xoTemEzVU0NHJcexZEsWKzNyje9xz2oIj4MY367pKnds\nxqoYqKqMyfalXwwYv5jWz4GCPRxU8WzPLiY82MyQbtHVDjt1SBeSYsLYkVPCor/3c8qQrk0a564v\n7qfHvrX8OvpZsjM0IUH7CbWYCAkyE2oxER4cRJ+ECIK8rdI98hrjl/vCmUa9uyWs3sMPlVRw59z1\nZBeV8+ZlR7ruSfhV2iJY/gzs/BVCOxrrMMZcW70w4OibYfhFsPQxo4Bg/Wcw4XYYewNYQv0/phYm\ngd1XJbnw/X2w7mMjFTHlSRh9OZgttY9VCvpNMv5kroYVz6GXP0PFurmEXL/UeL6f/JKWxQNfbcCu\n4Z9z1rI7t4Qbj++LqitnXpILPz1krCwcfLbX82/ZX8jUF5dTbvV83yBZXcmi4DtZ+cn/OPXWVxv9\nfSzdmsXm/YV0igrhjJHdan3dmZpwrsgkc7URFH24N1BWaaOw3EqQSREd1rw/8vG+bmrtlOR4x5GZ\nykrrGMCo3a8ZUIPMJq6Z0JsHv9nIK0szOHlwl7r/jb3YtH0XvTe+xRf2Y7j9107AWo/HnTUykWcu\n8JJ3NwfBqY/Du6cb5ZIT607N7cwp5op3/iAjy7gvc9NHa3jvyjHef3k0QMWvrxC8aCb2qERMJ/8X\nRl1a9+QlIg6mPA5jroFF/4af/gOpbxuVQQN9LCZopSQV44vSQ/D+mfDXZzD+FmPxzphrPAf1GraH\nDuDKkps5r/wBKNhLzrsXGzf+/GB3bgk3f/wndg3j+8ahFDz5w1bumrueijoCMT8/bJTXnfq416Co\ntWbWNxspt9pJigljdM8Yju4TxwkDOnHqkC6cNTKRMaPHsEyPYGTeQnYcbHzDKmf7gCuPSSYkyFzr\n68OSOhBmMZORVUx2djZkbfY5DZPrVurY2GDoqwZ1eASjVNMcDHtS+d3xS2usWxrG3XkpScSEW1i3\n+5Dnlbg+2vDDO4SoSn6LP5dzRydx+rCunDSoMxP6xTMmOZbhSR0INpv48s89rKxZieRJ8gQYco4x\nS87b4fGQ1TvzOOvlFWRkFTOgSxTxkSGsSM/hv99tbvT3UVNF2mLMi+7lB9toxhQ9yftMwRoU7v2J\n8f1g2sdw2bcQ2sGo0V/7kd/GFQiH1Yw9u6ic/NJKwoPNhAcHERFs9j5bKCswbrIc+Bsu/Aj6T/bp\nWoVllbz48zbe+nU7lTZNsHkQ91uv5IkDs9n3+V10Pf/pJn0vJRVWrnkvlUMllZwwoBOvX5rCj5sO\n8M85f/LZ6kz2HCrllYtH0yHM7ZfP3j+NGcnY66HzIK/XWPDXfn7LyKFjuIVvbzqmzq6I7+WfxQm7\nH+C97z+j1yVXNfh7WZ95iN8ycogMCeKio3p4PMZiNjG6ZwzLt2Wzbd0vxKON/LQPmrsPu7sGLVIC\nI8fdZRjsWcOqvJMAo37dk/DgIC4Z14vnf0rjtWUZ1fLwvjpQUEbfPd+whe7cPP18esR7Lm989set\nPPtjGrO+2ci8m4/x/v/kpIeNirCF98K06kFx/vp93PbpWiqsdo7rn8CLF41k8/5Cps1eyZvLtzMk\nMZqzRiY1+HtxZ83eTsXHl7DH3o277TeRV6p54OuNfLhqFw/+YzDj+vjwWiUfC1cvgo8vhK9vNH7h\nDj233qdorSmrtFNcYaW0wka51UZyfKSroilQ/BLYlVKnAM8BZuANrfX//HFef9Fa8/ovGTy2cEut\nG0/BQSYiHIG+Y7iF647rwz+GO1IBFcXw0flGQDz/PZ+Cut2u+XxNJo9/v4WswnKUgvNTkrjj5CN4\nZlESb6/ZyRV/v0nm4sEkHd/wIOj8fu6au57N+wtJjo/gmQtGYDYpTh7chU+vHceV76SyIj2Hc15Z\nwduXH0n32HCjBHPBnUaOv563y04lFVYemf83AHdMPqLeVrdHnTKd3NmPEZ/2Gfkll9Ih3Ps7GXev\nOXLrFx3Vg+jQup87JjmW5duyKdy20njAWS7ohXOJf3PfODWu4Wwr0ICKqKQU9Or3yCjKJyTIwtCk\nDnUeetm4nry2NJ2fNx9k64FC+neufwPvmub9tJirTGnMjbuOc+sI6mCUWM5dncnm/YV8uGoXlx3d\nq/4Td0iE4+40SmjTFkG/k9Ba8+rSDB5baMzKLzqqBw9NHUyQ2cSRvWJ5cOpgHvhqAzM//4t+naIY\nklj3910fe1khB18/m3CbndtNdzPnhpPIyCri/+ZvMn6BvL6SKUO7cO+UgSTFeJnBW8KMCdyH58EX\nM4zgPmgqYNwfeODrjWzYk09xuZWSChvFFdZaBUEDu0Yz+5LRxv+7AGlyKkYpZQZeAk4FBgHTlFLe\np4MtpKTCys0f/8mjCzZjs2uS4yPoEh1KVEgQJgUVVjt5JZXsOVTKxr0F3Pzxnzz/Uxq6osT4zb17\nlbFww4ec25+78jjrlRXcOXc9WYXljOrRka9vHM/j5w6nU1Qo/3fmENYMuINfbYNJWHo3ezb80qjv\nafayDOat30dEsJnZl1SflQ9L6shXNx5N/86RbDtYxFkv/8ra3YeMewOZf8BJDxlvN714ZUk6e/PL\nGNwtmmljPM+inY5IjCM16kROVKl8vmJDg76XnTnFfPfXPixmxRXje9V7rHMmG5691tg31Md7FVV9\nYpp/xu7KsfuaigFIHI2yltBfZTKqR4zHVJRTXGQI56d0BxreHKy0woZa9zFWbaLviVfWe2yoxcz9\npxn/jZ/6YYtv38/YG4zyze/uprK8lHu//MsV1O+dMoBHzhxSbeZ/8VE9uCClO+VWO9e+v7phr5mD\nttvY/OoldC7bzr/st/Lwlf/giC5RnDq0Kz/96zhuP6k/oRYTC/7az4lPLeWZRVsprfCyFiA4Ai76\nxEj1zb0StixkR3YxZ7+8gm/X7WV7djEHC8spKjeCekiQidiIYLrHhhEbEcymfQWc/sJylm0NXKWN\nP2bsY4BtWusMAKXUHOAM4G8/nLua1K2ZRFNM//5H+HT8juxirn1/NVsOFBIZEsRT5w/n5MFdXF/X\nWlNutbt++36/cT+PLtjES4s2MHnNDRxR/AfqrNdgSP03GXOKyvnfd5v5bLVRY9wpKoR7pgzgjOGJ\n1ao7zCbFUxem8M+3ZtEj8wZC517KgZjFdE7s5fNrsGxrlus/y9MXjKCfhxlbUkw4c68/mhs+WMPy\nbdlcM/tHloc/QEjSGBh2oddr7Mopcc2i/zN1sE9vK2PHX0HI91+T9dvHVB4/3Oee568ty8Cu4ewR\niXTtUH9FxfDuHQkOMtGnYgvlnU/E1zBd1Sem5WbsPqdiwHWvYLgpnS69J3k9/OoJyXy4aidfr93D\nbSf1J7Fj/a+b0xerdzJFL2NNSApHDuzv9fiTBxt591/Ssnnyhy389+xh9T8hKMRoOfHBOXzx8r18\nfGASIUEmnr1gBKcOrV3Fo5TioTMHs+VAIWt3H+LGj9bw/lVHNahf/qp372XsocU8ar+Yyy+9glE9\nYlxfC7WYueXEfpw7Oon/freZb9ft5bmf0pi7OpMHTh/EyYM7133PJSQKLp4L752B/ZNLeFzfSUbp\nYAZ2jeaJc4cRHxlCeIiZcEv1dG5+aSW3f7KWnzYf5LK3f+eOyUdww8Q+zX5vpyZ/3DxNBHa7/T3T\n8Zhfaa3JmXsrsR9O4oHnXuWrP/dQbq37N+/izQeZ+uJythwopE9CBF/dOL5aUAfjByvUYiYuMoTu\nseFcPaE3r180jFdDnmdA8e+81vFWDvU7q85r2O2aD1ft5ISnlvLZ6kyCzSaun9iHn++YyFkjk2qV\n7IGR+nnq8hN4Kn4WEbqY7LcuIDfftzr3nTnFrpult5zQt9b34y461MLbVxzJBSnduU1/iKUslwU9\nbvepXvmheX9TYbVz1shEUnr5NisefdRxZJh6MrnyJxb8tc+n5/y5K485v+/CpODaY71vRh1qMXNC\nt0q6qDwyQgb4dA1oucVJABHBZkKCTJRW2iip8PEmeWxvClQUI9S2aguT6tIzLoLTh3Wj0qZ54Kt6\ndpdyY7dr1i/9ii4qD9OI6T4FGqUUD/5jMEEmxZw/drM+85DX52R3mcAKy1j+cehDRoZn8/GMsR6D\nulNIkJlXLx5NQlQIKzNyeXSB7wuGFn/9NmN3vsoXtmMYfu59TOjnuc9Rt45hvDBtJJ9eO45BXaPZ\nc6iU6z5YzZXv/MGunHq6ZoZ2YMGIl9hs68bT9se5oUcmn103jiGJHejSIZToUEutew8dwiy8fmkK\nt07qh9bwxPdbuP6DNRSV+6dgwlf+COyefkJq/aQppWYopVKVUqlZWQ1/i1JutZPW5woKieTB3HtY\nO/d/HP3oTzy+cHO1lqZ2u+b5n9K48t0/KCizMnlQZ766cTx9O3mv18ZWyYl/38Pxag2Pma7hf/uP\n5OxXVrAzp/ay+b8y8znrlRXc9+UG8ksrmdAvnu9vO5a7TxlAZEj9b4TCg4OYdfX5PBVxG4Ntm0l9\n+UqKyurvcV1SYeXa91eTX1rJiQM6cesk7zMui9nE/wZt56Kgn5ltO50bfrbzypL6tztbsuUgP246\nQESwmXtO9T14KpOJQ/3PY6RpGwsXL/UabMqtNu6aux67hmuO7e3xnYcnp8UYLYNXlPXyeWzZLbQ4\nCYxgGN/AWXt+qZU1tt6MNKdXm3HW577TBhIVGsTPmw/yzbq9Xo9fvOUgxxR/Tz6RjDjxAp+uAdC3\nUyRXjO+F1jDrm43Y61kctedQKee/+hu3F06nUgUzJ/4NRnXz3numS4dQXpk+CotZ8favO/jCh9W1\nPy9bwpFrZrLW3hvrlGc4bXjtEtmaxiTH8u3Nx/DwGYOJCg1i8ZYsTnpmKc//lFarVYPWmhd/TuOG\nL3YwvXwm+WFJ3Jk3i8j9v3u9jsmkuHVSf968LIWo0CAWbtzPmS/9SnpWyzVM80dgzwS6u/09Caj1\nk6a1nq21TtFapyQk1NNBsA6hFjM3XXAanW7/lf2dj2WW5T3urXyeN5ds4tjHF3P1u3/w498HmPH+\nap5etBWAOyb359WLRxNVzw05l8oy+PJao2/GKf/jklseZkCXKDKyijnr5RWs3mmUl+WXVvLvrzdw\nxkvLWbf7EF2iQ3l5+ijeu3IMyfXcjKqpY3gw1153G+8Encfk8h+Y+8q/KS63UlxuJbe4gn35pezI\nLmbrgUL+ysznjs/WsXl/Ib0TInjmwhEe3w3Ucmg36pubodsooqbMMtrgLtzMYws3ewy8FVY7D31r\nZNBuObEfnaIbtlBj0MlXY8XEsJwF/LEjr95jX1qcTtrBIpLjI7jNh19STiPNGZTrIL496HtFSEv0\nYncXF9mwRUp/7Mhlrb0P/VQmoXbf+q53jg7l/tMGAkbA9XatOcvWM9m0msyk0wgK8S1143TLif2I\njwxhza5DfPmn5178GVlFnPfKCjKyi4npmgxTXyDk4Hr4+SGfrpHSK5ZZUwcDMPOLv/hh437+3ltA\nelYRu3NLOFhYRn5pJWWVNn5Zt5U+P82ghFDWj3+Z88f5/vNjNikuGdeLn/81kbNHJlJutfP0oq2c\n8uwyV068wmrnzrnrefKHrSgFN58+lk43fo+KTjRuqu5c4dO1ThzYmW9uOsZ1v+vMF39l0d8HvD/R\nD/yRY/8D6KeUSgb2ABcCF/nhvB5FdIgl4rov0cse55wl/2VM+AGmF93Mj5vgx00HAYgODeK5aSM5\n/ohOXs7msGsVfHMTZG+FSf+BsdfTDZh7/dHc9NEalmzJYtrrq7hyfDJzV+8mu6gCs0lxzTG9+Oek\n/l5n6HXpFB3K8dc9w9KXdnDxoVe57D8d+NU+tM7jI0OCmH1JSr2VIy42q3FX326Hc99kemxvIsPD\nuP3TdbyyJJ2C0koePmNItV8Qb/+6nYzsYnrHR3DFeO+7EtUUGtOVbbHHcHbOLzy4LI0xyUd5PG7T\nvgJeXmw0jvrf2UMJtdR9s7CmrkUb+Uv3ZP2+UgrKKn16LZx9W1qi3NG4TsNKHldtzyHN3hcT2ljZ\nnDzBp+edn9Kdb9ftY/m2bP7z7d+8MG2kx+M27s2n8675hFgq6Tlphm/fhJuoUAszTx3AHZ+t47/f\nbWby4M7VJksb9+Zz2Vu/k11UweieMbx1+ZHGDf19y2HFC9B7otFB1IvpR/Vkw558Pv59NzPer92n\nHiAIK29ZnqCLKYdPBr3CpSePa/D3A5AQFcLTF4zgvJTu/PvrDaQdLOLSt35nytAu5BVX8ltGDmEW\nM89dOILJzpTnZd/AO6cZf8bdaGw+E1x/5UtyfARf3jCeO+euY8Ff+7nmvVRuP6k/t5zYr1Hj9lWT\nZ+xaaytwE/A9sAn4VGu9sannrZfJhJo4E6bNoTv7WRI9i2fHFNA9NoxhSR349uZjfAvq5UWw4C54\n62Sjy+HFnxv9TxwiQ4J449IULh7bgwqrnVeXppNdVEFKzxjm33IM9502qNFB3alnfBRdLn+XXaYk\nZlueZlxwOjHhFrpEh9IjNpx+nSIZkhjNuN5xzL50tG8pJTDaqe5aYWxvFmvkr88YkchrF48mOMjE\nh6t2cduna119ww8UlPH8T8Y2bf/+xyCCgxr3o5FwzBV0Voco2/ojOzxsQG212bn78/VY7ZpLxvZs\nWC223YZ531r2RgzCriF1h2+LdFx17C02Y3ekYnxpBIbRH2ad3XGPwdkDxwdKKf579lDCLGa+Xbe3\nztngm8u3c655GQfC+hDZ07dFXTWdPTKRkT06kl1U7vo5AePf4MLZK8kuqmBCv3jev2pMVZXW5P+D\nToPgy+uMFtE+mDV1MJeO68nw7h0Z0CWK3vERJHYMIz4yhI6hiueDX+ZY818s7HU3l5x/XqO+F3fj\n+sQx/5YJzDx1AGEWs2vtRkJUCJ9eO64qqIOxq9SMJUbPmRUvwKvjYbv3yraIkCBeumgUM08dgElB\nTAPLgRtFa93if0aPHq39Jmur1i+kaD0rRusVL2ltt/v2vLQftX56iNYPdtB6/h1alxXUeajdbtdv\nL8/QU55bpj/9Y5e22Xy8RgPY8/dq/dxIrR9N0nrPmqadbMevWs/qqPXnMzx+ecW2bD3oge90z7vn\n6Svf/l2XVlj1rXP+1D3vnqeveuePpl27slwXPdRdf3v/ZP3vr/6q9eXXlm7TPe+ep8c9+qMuKK1o\n2Ln3b9D6wWj9zXtP6553z9OPLvjb61Psdrvud+8C3fPuebqk3Nqw6zXSowv+1j3vnqdf/DnN67GF\nZZU6eeY83fue+dr2zHCt50xv8PXe/CVD97x7nh7zyCJ9qKT6a7o/v1Sfcu9rWj8YrXMXPdXgc7tb\nv/uQ7jVznu5zz3yddqBAL91yUA+43/g5uu79VF1W6eH1PfC31g930vq9s7S22Rp/cZvN+Hl+MFrb\nf32+8eepR2Zeib7pozV6+usrdWZeSf0Hpy/R+tlhWj8YrfU3/9S69JBP19i8r0DbfY1RHgCp2ocY\n2/ZbCsT3g6t/giNOhe/vgddPMBbirHkf9q2v3Tq3NA++usFo2xkUYmz3NeUJo7ypDkopLh+fzPxb\nJnBeSnff8tsNpKK7Gm/1wjrC+2fBgUa+6SnNg8+vMZpjneZ5E4RxfeL46JqxdAy38NPmg5z98gq+\n/HMPwUEm/n16E5cgBAVTPvAcTjKlsjB1M/klVTeFt2cX89QPxv2PR84e6tu9D3eZxmy208DxQNXe\noPUpLLdSYbMTEWwmLNj3lE9TOBuN+ZKKSd2Ri13D0MQOmLqnVGvh66vLju7FqB4dOVBQzn9rVJW8\n99sOzlDLsGEiZuzFDT63u6FJHbjwyO5Y7ZrrP1jDVe/+QWmljfNGJ/HCtJGe6+87DYSTHzU6QK58\nqXEX1hrm32Y0Szv+ftTRjWgx7INER/XMB1cf5b2EtPdxcP1vMO4mo6neS2Nh6/e1jys8YCzYWvYE\nfHIJR3xyDMrDdoj+1vYDO0BoNJz/Ppz8X2Ol2J8fGjnz1ybAo93g1QnGEuFlTxgbNq+bA8fcDtct\nh56Ny9E1iw5JcOk3EBRm7EPqYQf7emlt9J0u2m/sT1rPL6vh3Tvy6bXj6Bwdwt/7jB4v1x7bmx5x\nTV8tFzv+ckKUlZPsy/n4j12AUa008/P1lFvtnD0y0ff7H+72pEJoRwYNHoFJwV+OFYD1yWmBnZNq\nakhP9lVu2+CRONrYlajAe5WLO7NJ8fi5wwg2m5jzx25+3ZYNGAuS5qzczlnmXyjsfrzRZreJ7ph8\nBNGhQaQdLKLSprnqmGQeO2dY/S0HUq40duj68T8N33VLa1h4D6x+Byb8y1jd2loEh8PJj8BVPxoT\nso/ONxY0/fSwsZnHk/3hqf7w4bnw8//B/r+MFt+m5u/k0j4COxj12eNugKu+h3sy4abVRnAbd4PR\n3nXLd8aLG9kZZiyGSQ+2zvacscnGzB3g3amQ24BNjFe/A5u+gRMf9Gm5ff/OUcy97mgGdIlicLdo\nrp/Yx+tzfNJlGEUdB3CueRnv/LqDSpudj37fxartucRHBvNAY98V7FkDiaOJCgtmSGIHbHbNml31\nV9+0VB92dw1ZpOTa37R3bFXvm0bM6Pp2iuKWE40NU2Z+sZ6SCiufr8lkaPkaOqtDdBh3WYPP6Ulc\nZAizpg4mMiSIf53Un/tPG+j9HaxSMPUF4xfL3Ct9359Aa6MT6apXjFWtJzzQ9G+gOSSNhhlLYeI9\n8Pc3RjO0gr1GG+OT/wuXL4CZu+Cfa43WJA3ZqaqR2mcTMJPJ2BUovm9VEx+toTjbWIZuapm35I0W\n3w8u/dq4+/7uVLjyO2M2X5+Dm6p6Yo+7yedLdY8N57t/TsCu8V/jIqUIH3MJI364j8jCbbzxSy9e\nclTB/GfqkHr7ztSpvMjYNWnAaYAxw12fmc+qjNw6F6aA216nLVQRY1zLt3LH0gob6zMPoZRR7oe5\nA5gsRspp4D8afN1rj+vD/L/2s2lfAU98v4WlW7K43byMiuCOBPc/tVHfiydnj0rijBGJDft5CY81\nNqF593SjYOGsV7w/Z9mTsPxpGH2Fkc5p4dWbDRIUbPRgOuo6I8XrpS99c2s/M3ZvlDIa7bf2oO7U\neTBc8hWU5cO7/4DC/bWPsZYbbVJ3/gZzrzJSL2e+2uDdcJRSfu9GZxp2AXYVxLnmZTy2cDNF5VZO\nHtyZKUPrXi1br31rQdtdy+/HuPZBrb+trDMd0lI17Ma1nFUx9c/Y1+zKo9KmGdQ12ijbtIRClyGN\nmrGDsSDtiXOHYTYZC31ysg8w2ZxK0IgLfNrwuyEa9fPSazwceyes+wh+fd5ITZTket5Wb8ULsPj/\nYPg0OO3p1h3U3YV1DHhQh/Y6Y28vuo1w9Ks40wjuPcdD4T7jbV7BXijJrjpWmeCiTyGqc+DG6y4y\nAd1vMudsXc4T1guICA3h4TOGNL5nhjPYOQN7r1iUgnW78ymrtNVZC9/SpY6Aa1eg3OIK7HZdZ6qi\nKr/uVvKZmGI0bLPbGjUJGZLYgWsm9ObVpen8w/wbwVhh5PSGfxPN5di7YMdyWPQALHI8FhQKUV2N\ncsLobsZ9srUfwuCzYOqLjdu27zAngb216z7G6DT32WVG/jzK8cOfOKrq8+iuxjZ9HevvwtjSzCOn\nk7B1Acea13PGGVc0eCVrNZmpRqWPY2vBDuEWBnSJZtO+Av7cdajOfttVDcBaLhUTHGQiOjSIgjIr\n+aWVdaaeVjk2sTiqt1t/mKQU+ON1yNriU898T26d1I+lW7OYdugXbPGDMXfx0ryrJZmDjHei+9ZW\nTVCcN4wL9hmbrBcdMHb3Onu2cbxoMHnV2oLkCXBnett5O+rUbzI6PI7ZPbZiGdnEvnB7VkOPsdUe\nOio5lk37Cli1PafOwJ5d3PIzdjDSMQVlVnKKyz0G9rJKG3/uNppqHeneaM25K9Se1EYH9lCLmW/O\nj8Xy2jYY1Qpz00HBxoSlLlq3vjG3MfIep61oiz/oQcGooedjSfsOin3YYq0uBfugYE+tHZOc/dkX\nbthPYR1N1Jwz9pboxe6uql+M5zz7vPX7qLDaOaJzVPUNnWP7GP3yM31fgeqJZf1HRlnd0PObdJ6A\naIs/662MBHbRvEZfDrZKWPFc48/hzK8nVQ/s4/vFEx8Zwub9hZz7ym/Vunw6BSLHDlWpn5olj1pr\nXl+WwZ1z1wFwZs13MiaTMWvfs6bxFy/KMrZAHDjVKBgQhx0J7KJ5dRoAw86HVbONmXdj7Ek1Zp9d\nqjdIiw618Pn14+iTEMGWA4Wc+dIK1u2u3jPc1Yu9BXPs4HmRUqXNzn1fbeCRBZvQGu48+QiuO85D\nH/rE0XBwo7E1Y2MsfxqspXD8vY17vmjzJLCL5jdxJtgrjcZkjZGZCp2HeCwj6xkXwRfXj+foPnFk\nF5Vzwezf+M6x0YfVZievpALVUo2X3Lj2PnXM2AvKKrnynT/4aNUugoNMvHjRSG48vq/nKqHEFKO0\nc+/ahl84PxP+eANGXGSshxCHJQnsovnF9oZRlxorYxuykhaMsr+9a2ulYdx1CLfw7pVjuCClO2WV\ndq7/cA2vLEknt6QCrSEmPLj+Je/NwH3v0925JZzz8gp+ScsmLiKYOTPGcvqwejaGcN1AbUQ9+9LH\njY/H3d3w54p2QwK7aBnH3mmkU5Y+1rDnZW+FisKqYFcHi9nE/84ZykzHrk+PLdzM7Z8YeeyW2Ou0\nJmfqZ/XOPM586VfSDhbRr1MkX9043vsuSZEJRulqA1r4ApCTDn9+YPRmaWWlr6JlSWAXLSO6G4y5\nxmjAdnCz789zVock1j1jd1JKcd1xfXj14lGEWkwsdzTDaukbp+7X3Ly/kJxio1f55zccTfdYH5us\nJTai0+PiR43l7BP+1cDRivZGArtoOeNvg+BIWPyI78/ZvRJCOkBcX5+fcsqQrnwyYxwJUcasOSGq\n5Zu9Oa8NcNFRPXjr8s9LuKYAAAdGSURBVCN92/nKKSkFCjLh0G7vx4KxPH/DXBh7vV+6OIq2TQK7\naDkRcXD0TcYKWl/K+da8b7RgPuLUBi8rH969I1/dOJ5Lx/VkxgQPlSfNrHd8BLec0JfHzhnKI2cO\nwdLQHH/fSWAOgTkXGf1UvPn5EaP+vZl6lYu2RWkvu8k3h5SUFJ2a2rQFGKKNKiuA54ZDt5FwyRd1\nH7f6Xfj2FuhzIlz4Uetssdzc0n40AntCf6NPf3is5+N2/w5vngQn/lvSMO2cUmq11tprXlJm7KJl\nhUbDMbcZO+rsWO75mNXvGEG976TDN6gD9JsE0z6CrK1G++a6Vu/+9BBEJBgtY4VAArsIhDHXGN38\nfnq4dsvW1Lfh239Cv8lwwYeHb1B36jsJpn0MOWnwnofgnrEEdvwCE+6A4IiADFG0PhLYRcuzhBnl\nj7tXGvtBOv3xJsy7FfqdDBd8IEHdqe+JMG0O5Gwz2jcXO9o1O3cYik6ClCsCO0bRqkhgF4Ex8hKj\nDe/PD4HdbqyWnH879D8FLnjfKNsTVfocb7Rvzs0wgntRFmxZYCximjhTXi9RTZMCu1LqPKXURqWU\nXSnlvdBYCKegYJh4r1GmN/cKmP8v6H+qsSekBCnPek90BPftRnD/8T9GGejwaYEemWhlmjpj3wCc\nDSzzw1jE4WbouZAwEP7+Co6YAue/K0Hdm97HwfRPjS0Rs7fA8ffJZhSilib9RGitNwGN3+5MHN5M\nZjjzZdjynZFz9/O+nO1W8rHGZufpP8OgMwM9GtEKtdiveqXUDGAGQI8e0sdCOCSOMv6IhulxlPFH\nCA+8Bnal1I+Ap63l79Naf+3rhbTWs4HZYCxQ8nmEQgghGsRrYNdaT2qJgQghhPAPKXcUQoh2pqnl\njmcppTKBccB8pdT3/hmWEEKIxmpqVcyXwJd+GosQQgg/kFSMEEK0MxLYhRCinZHALoQQ7UxANtpQ\nSmUBOxv59Hgg24/D8ScZW+PI2BpHxtY4bXlsPbXWCd5OEpDA3hRKqVRfdhAJBBlb48jYGkfG1jiH\nw9gkFSOEEO2MBHYhhGhn2mJgnx3oAdRDxtY4MrbGkbE1TrsfW5vLsQshhKhfW5yxCyGEqEebCuxK\nqVOUUluUUtuUUjMDPR53SqkdSqm/lFJrlVKpAR7LW0qpg0qpDW6PxSqlFiml0hwfY1rR2GYppfY4\nXru1SqkpARpbd6XUYqXUJseWj/90PB7w166esQX8tVNKhSqlfldKrXOM7T+Ox5OVUqscr9snSqkW\n30mlnrG9o5Ta7va6jWjpsbmN0ayU+lMpNc/x96a/blrrNvEHMAPpQG8gGFgHDAr0uNzGtwOID/Q4\nHGM5FhgFbHB77HFgpuPzmcBjrWhss4A7WsHr1hUY5fg8CtgKDGoNr109Ywv4awcoINLxuQVYBYwF\nPgUudDz+KnB9KxrbO8C5gf6Zc4zrduAjYJ7j701+3drSjH0MsE1rnaG1rgDmAGcEeEytktZ6GZBb\n4+EzgHcdn78LBGRPtTrG1iporfdprdc4Pi8ENgGJtILXrp6xBZw2FDn+anH80cAJwFzH44F63eoa\nW6uglEoCTgPecPxd4YfXrS0F9kRgt9vfM2klP9gOGvhBKbXasQ1ga9NZa70PjCABdArweGq6SSm1\n3pGqCUiayJ1SqhcwEmOG16peuxpjg1bw2jnSCWuBg8AijHfXh7TWVschAfv/WnNsWmvn6/aI43V7\nRikVqF3UnwXuAuyOv8fhh9etLQV2Tztmt5rfvMB4rfUo4FTgRqXUsYEeUBvyCtAHGAHsA54K5GCU\nUpHA58CtWuuCQI6lJg9jaxWvndbaprUeASRhvLse6Omwlh2V46I1xqaUGgLcAwwAjgRigbtbelxK\nqdOBg1rr1e4Pezi0wa9bWwrsmUB3t78nAXsDNJZatNZ7HR8PYvSoHxPYEdVyQCnVFcDx8WCAx+Oi\ntT7g+M9nB14ngK+dUsqCETg/1Fp/4Xi4Vbx2nsbWml47x3gOAUsw8tgdlVLOPR8C/v/VbWynOFJb\nWmtdDrxNYF638cBUpdQOjNTyCRgz+Ca/bm0psP8B9HPcMQ4GLgS+CfCYAFBKRSilopyfA5OBDfU/\nq8V9A1zm+PwywOeNyJubM2g6nEWAXjtHfvNNYJPW+mm3LwX8tatrbK3htVNKJSilOjo+DwMmYdwD\nWAyc6zgsUK+bp7FtdvtFrTBy2C3+ummt79FaJ2mte2HEs5+11tPxx+sW6DvCDbx7PAWjGiAduC/Q\n43EbV2+MKp11wMZAjw34GONteSXGO52rMHJ3PwFpjo+xrWhs7wN/AesxgmjXAI3tGIy3veuBtY4/\nU1rDa1fP2AL+2gHDgD8dY9gA/NvxeG/gd2Ab8BkQ0orG9rPjddsAfICjciZQf4CJVFXFNPl1k5Wn\nQgjRzrSlVIwQQggfSGAXQoh2RgK7EEK0MxLYhRCinZHALoQQ7YwEdiGEaGcksAshRDsjgV0IIdqZ\n/wfIUJFhjYQYHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\zhaon\\Anaconda3\\lib\\site-packages\\numpy\\core\\numeric.py:501: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd8m9W5x79Hy7Isz9iOs+OE7AkJ\nSSAFwiXMC4Q9ymW1gVJa2kJ7gQ4KnfQyuggr7NmySgshhJlAgBAyyB5kOMOJ99SwrHXuH68kL8l2\nbMfSK53v5+OP7FfvOLbf99HvPOcZQkqJQqFQKFIHQ7wHoFAoFIr+RRl+hUKhSDGU4VcoFIoUQxl+\nhUKhSDGU4VcoFIoUQxl+hUKhSDGU4VcoFIoUQxl+hUKhSDGU4VcoFIoUwxTvAUQjPz9fjhw5Mt7D\nUCgUCt2wbt26aillQXf2TUjDP3LkSNauXRvvYSgUCoVuEELs7+6+ytWjUCgUKYYy/AqFQpFiKMOv\nUCgUKUZC+vgVCkXq4fP5KC0txePxxHsoCY3VamXo0KGYzeYen0MZfoVCkRCUlpaSmZnJyJEjEULE\nezgJiZSSmpoaSktLKS4u7vF5unT1CCGeFkJUCiG2tNr2ihBiQ+hrnxBiQ4xj9wkhNof2U2E6CoUi\nJh6PhwEDBiij3wlCCAYMGNDrWVF3FP+zwCLg+fAGKeXlrQbyINDQyfGnSimrezpAhUKROiij3zV9\n8TfqUvFLKT8FamMMQACXAf/o9UiOMnurnLz45X48vkC/XtfrD/Ly6gPsLHf063WTlqqdsPpxOPx1\nvEfShmpnM8+v2keVo7nvT35wjfY71+zp8JaUkqWby/h8t9JWiu7TWx//SUCFlHJXjPcl8L4QQgKP\nSykX9/J6PcLrD/I/T67mcIOHbyoc/HbB5H679iMrdvPXD3eRYzPz2R3/hT1NLav0mG/eg39eBUEf\nIOCc+2HWDfEeFQC3vrKBlbuqeWdTGa9874S+O/Fnf4EP79G+N1nhqteg+OTI20s3l/ODl9cD8MGt\nJzNmYGbfXTsF+fvf/86jjz7K9OnTqaqqorq6mp///OeUlZVx4403YrPZ4j3EPqG34ZxX0rnanyul\nPA44G/iBEOLkWDsKIW4UQqwVQqytqqrq5bDa8vnuag43eMjLsPDGutJ+U/1ef5AXVu3HZjFS7/ax\nYmdlv1w3KWmqg3/dCIUT4Adfwbiz4d3b4cDqeI+MKkczK3dVk2Exsrqktu9md3s+1oz+5Ivh5tWQ\nMwLeuAG8rsguL63eT6bVhNEgeH19ad9cN4V55JFHWLp0KT/+8Y/x+Xxs2LCByy+/nL/+9a+43e54\nD6/P6LHhF0KYgIuAV2LtI6U8HHqtBN4EZnWy72Ip5Uwp5cyCgm6Vm+g2q0tqMRsF918yFZc30G8G\neMPBempcXu6/ZBr2NBNf7Knpl+smJaseBk8DXPAoFIyDi56AzEGw9GcQ7F/3XXs+260JlUf+ZwZC\nwNLNZb0/qd8LS/8X8kbDgkegcDwsWATOcvhKmzhXO5tZtbeG78wt5rjhOawpieqRVXSTm266ib17\n93LOOecwd+5cNmzYwPTp0/nb3/7G4cOHOfXUUzn11FPjPcw+oTd+h/nADillVJkhhMgADFJKR+j7\nM4Df9uJ6PearkhqmDMnm5LEFpJuNfLm3lrMmD+qX6wKcOHoAs4rz+Eo9mD0j4If1L8CYM6Ao5KZL\ns8P838C/FsKOJTBxQdyGt2ZfHVlWEycdk8+kwVms2lvDrb096eZXoWY3XPkKmK3atmGzYORJsOZp\nOPFHfFVSi5RwyrgCPL4Az3y+D48vgNVs7O3V485v3t7KtsONfXrOiYOzuPu8STHff+yxx1i2bBkr\nV65ky5YtPPDAAyxZsgSAv/zlLyxfvpz8/Pw+HVO86E445z+AVcA4IUSpEOK7obeuoJ2bRwgxWAix\nNPTjQOAzIcRG4CvgHSnlsr4beveQUrKj3MHUoTmYjQaOHZ7TbwZ4dUkt44syyc2wMGFQJiXVLrz+\nYL9cO6nY85GmdI+7pu32yRdB7kj4YlFchhVmd4WTsQMzMRgEJ4wawIYD9TT7ezELkVKb4QycDGPP\nbPvejOug4QAc/IqvSmpJNxuZMiSb6cNy8AaCfFOhgggUXdOl4pdSXhlj+3VRth0Gzgl9vxeY1svx\n9ZqKxmbc3gCjCjIAmDkil0XLd+P2+rFZjt5Cqz8QZN3+Oi6dMRSAYwrtBIKS/TUutQB3pGx/C9Ky\nOxpBgxHm3Kz5+g9+pSnifkZKyTeVDs4OzSCnD8vliUAJuyqcTB6S3bOT7vkYKrdpbq32oXtjzgCD\nGXYsYVPpOUwZmo3ZaGB0oR2AkmoXU4fm9OZXSgg6U+aK3pP0tXr2VjkBGJWvPRiThmQTlBzRApwv\nEGTd/tojWhTeV+PG7Q1EHsJjCjRjv7vS2e1zKNDU7+6PYPSpYIySoj79KrDYYf1zR3ccwSCUb4by\nLdr3IaqdXurdPsaEDO/kIVkAbDnUWWpLW6SUfH2gjjqXV9uw4WVIz9UWddtjzYJRpyB3vMPO8kYm\nDtKuNzxPizYpqXZ1PEbRazIzM3E4kmc2lfyGP/QghBV/+EHZVtY9/2EgKPnOs2u4+NFVXPLYF912\n1WwPnX9C6Hrh6+9VD+aRUbEVHGUw5vTo76fZNf/+1v+A9yhFXZSug0dmw2Pfgsfmal8VW7W36rRr\njhigGd5huTYy00xsOdx9w3/PW1u58JEvOOOvn1JVUwM7l8KkC8GUFv2AMWci6krI9VUwvkgTFFaz\nkcHZVvap++uocOONN3L22WcnzeJu0hv+0romzEZBUZa2QDY0N51Mq6nbC0dvbTzEyl3VnD25iC2H\nGnlpdfd6Hewob8RkEIwu1Ax+RpqJ7HQzZQ1NPftFUpV9K7XXUZ08cNOuBK9DW+Ttaw6shufOBV+T\n5no5fxG4a+CpM6B8M+UNWup8UbZ2fxkMgomDs9jazftry6EGnlu1n/kTCmlw+/jwzafB54Ypl8Y+\naORcAGaJ7YwPCQuAkfkZ7KtJnpDDeLBv3z7y8/OZN29eZGEX4JZbbmHHjh0sX748jqPrO5Le8Fc2\neijMtGIwaL5SIQQTB2V1W/E/+8V+xg608/C3j2PGiFxe+HI/Usouj9te5mB0gZ00U0uExaBsK2X1\nqvLgEVG6FrKGQvaQ2PuMmAvZw2HDS317bXctvHoNZBbBwo9g+rfhuKvhhuVgzYaXr6C6RsuYHZSd\nHjls0uBstpc1Egh2fZ+8tHo/6WYjf758OpfMHMqgA+8QzBoKw+bEPqhgAh5TFrONOxg70B7ZPCg7\nPfJBpFB0RtIb/gqHh8KstlPmiYOz2FHmINjFg7m/xsXGg/VcMmMoBoPg0hlD2VvlYsuhrj80dpQ1\nMn5Q20XcwTnpHFYP5pFRugaGzuh8H4MBpl0Oez+Bxj6IoQ+z5FZN3V/6HGQObNmePQQuex4ch5m0\n9UEsJgO5tpb1h/GDMvH4ghyo7Vx9B4KSdzaVcfbkIrKsZq6cZONbYiO7Cs/UfqdYGAx8Y5nEiaZv\n2gQoDMq2UuVs7tYHjiK1SX7D39jMwExrm23jizJp8gW6fDA/C9U/OX1iEQDzJ2oP/6e7Os8sbnD7\nONzgYXxRVpvtg7KtytVzJDiroH4/DJnZ9b5TLgMkbP1X31y7ZCVs+zeccjsMmtrx/aEzYc7NHFf1\nJqfaD7QpnDUh9H/f0cWsckd5I40ePyeP1RIWJ9UvxySCvOHruuTDV/7RDJOHoak+sm1gtpVAUFLt\nPAr1ghRJRQoYfk/E/xpmXPjB7CKy56uSWgoz0xgZWrjLt6cxeUgWn3zTueHfXq498OGFtzCDc9Kp\nd/to8sY301Q3HApV8h56fNf7FoyFQdNg06u9v24wCO//CrKHwYm3xN5v3p00GLL5QbBt1ZIxA+0Y\nBGzv4v5avVfLJ5lVnAeAYfNrlKcV88r+LPyB2EEEbq+flc6Q66t8U2R7eB1LuXsUXZHUht/t9ePw\n+Du4esYOtCNE5yGdUkpW761lVnFeGzV30pgC1u+vw9nsj3ls+LztXT0Fmdo4lCLrJqVrQRg1g94d\nplwGZRugOlbNwG7yzbvaeU79JZjTY++XlsmzxouZ6v1aczOFsJqNjMzPYGd554p/dUkNw/LSGZyT\nDnX74eCX1I2+gAaPn82dhIN+U+FkS3Ck9kPZxsj2iOFvVIZf0TlJbfgrGzUD297VY7OYGJFnY2dF\n7AeztK6J8kYPs0NqLMy3jsnHH5Ss3Rc7+3dnhYMsqynyIIYZkGEBoCYcr63onIqtWl0eSzcrIk6+\nCBCw+bWeX1NKrSJmzvDOI2uAYFCy2H0KjZZCWPGnNu9NKMrqdEYppeSrklpmjRygbdjyBgD5J2j5\nkuv218U8dkdZIzVk47cPbmP4B2ZrwqLyaJSGViQVSW34w8qnvasHYFxRJjvKYj+Yq0NlHY5vZ/iP\nHZ6DySAi70djZ7mD8UVZHRomDLBrD2atSz2Y3aJyGxSM7/7+WYOh+CTN8Hcj8ioq+7/QFpRP/BEY\nO8/srm/y4QqY2T7yGjjwhVY3P8T4okz217hxxZgZ7qp0Uuf2MXtU6P7a/DoMm03BsHEMy0tn7b5O\nDH+5gwyLEePgaVDW4urJtWnCotaphEU8efbZZ/nhD3/Y5T6HDx+O/Lxw4UK2bdt2tIcWIakNf0XI\n8A/M6pgIM64oi301rpjZuGtKaslONzO2sK27xmYxMXVodsx6P1JKvil3MK6oY1mGsOKvVg9m1zQ7\ntYXdwolHdtyUy6B2Lxxa37Prfv5XsOXDsf/T5a61oZlb1djLwZoDX/wt8l74/x+rds7qvVoBvznF\nA7SZTeXWyAxj5og81u6vjRk2vL2skXFFmYiBE6F2j1bJEzAbDWRZTdS51f2V6LQ3/E8++SQTJx7h\nvd4Lktrw14QM7ICMjoZ/fFEmQQm7KqKXUPhqXy3Hj8yLxP+3ZvaoAWwqrY+6SHuovglHsz+64beH\nFJly9XRN1U7ttXDCkR038XwwpvXM3VO+BXa9D3Nu6ty3HyJsYLOzc+H478L2JZEuWeGM7VjuntUl\ntRRlWRmWl64tSBtMMPECAGaOzKXa6WV/lGQsKSXbyxq1xK38cRD0ax90IfIyLMqV2AsuuOACZsyY\nwaRJk1i8WCt/bbfb+eUvf8m0adOYM2cOFRUVALz99tvMnj2bY489lvnz50e2h3E4HBQXF+Pz+QBo\nbGxk5MiRvPbaa6xdu5arrrqK6dOn09TUxLx581i7VgtmWLZsGccddxzTpk3jtNNOOyq/Z1K3g6p3\nexECstI71ngJR9zsKG9kytC2xbQqGz2UVLv49qzhUc87qziPR1fsYf2BOuYe07ZMa2RhN4rht1lM\nWM0GatTibtdUhqa9R2r4rdkw9gzNZ37G77t017Th879pdX+OX9it3cMf4Lk2C8z6nlYldNUiOPcv\nDMlJx55mihrSKaVkdUktJ44egJBBzfAfczrYtbDOmSM098/a/XWMzM9oc2xZg4dGj58JRZna+gdA\n9U6tXj+a4a9LBsP/7p1abaS+pGgKnP2nTnd5+umnycvLo6mpieOPP56LL74Yl8vFnDlz+MMf/sDt\nt9/OE088wa9+9Su+9a1v8eWXXyKE4Mknn+S+++7jwQcfjJwrMzOTefPm8c4773DBBRfwz3/+k4sv\nvphLL72Uhx9+mAceeICZM9uGKldVVXHDDTfw6aefUlxcTG3t0akknNSKv87tIzvdjDGKah8xIAOr\n2RBVkX21r22YXXtmjsjFIIjq5w+fb2wUww/a7EMpsm5QuV1rNZg78siPnXIZuCph36fdP6Zuv/Zh\nMeM6rUBadw4J/R/zMixagte0K+Drl8BZhcEgGDvQHvX+Kql2UeVoZnbxACj5FByHtQS0EGMK7WRZ\nTVEDCHaUt6oBlT9G2xieHYXGomaUPefvf/97RNkfPHiQXbt2YbFYOPfccwGYMWMG+/btA6C0tJQz\nzzyTKVOmcP/997N169YO51u4cCHPPPMMAM888wzXX399p9f/8ssvOfnkkykuLgYgLy+6DeotSa34\n69zeyIJXe4wGwZjCzKghnV+V1GKzGJk0OCvKkZBpNTNpcHbET9uaneUOhuSkk2WNUkkSzd1To3z8\nXVO1XVO0hh40FRlzhlbGedNrMPq/unfMF38HYYATftDty9S6Wyl+0GL+1z8Pqx+F037N+EFZvLOp\nDCllm4X+sGCYVZwHn/0hVHL67Mj7BoNg5si8iABpzfayVsLCYtaij1oZ/lybpVuZ5QlPF8r8aLBi\nxQo+/PBDVq1ahc1mY968eXg8Hsxmc+T/ZzQa8fu1BftbbrmF2267jfPPP58VK1Zwzz33dDjn3Llz\n2bdvH5988gmBQIDJkzvv993+XjlaJLXir3f7yLFFN8CguWOiKv6SWmaMyMVkjP3nmV2cx9cH6zss\nDm853MCEQbHr7Q/IsFCjonq6pnoX5I/t2bFmK0w8D7a/3b2KnY4KrcPX9Cu1yKBuUufykm42km4J\nfTjlj9Eqha5eDE11TCjKpKHJ1yGu/quSWvLtFkbbvbDtPzD5wpYuWyFmF+ext8pFpaPtsdvLGhma\n20pY5I9rq/jtFmrd3m7Vk1K0paGhgdzcXGw2Gzt27ODLL7/scv8hQ7REuueei10W/JprruHKK69s\no/ZjlXk+4YQT+OSTTygpKQFQrp6e0JniB61mT7WzuU0ZhVqXl50VDmaN7HyKNas4D68/yKbSlkSb\nereXvVUujh0e21WQl5Gmwu26wueBhlKt32xPmfZtrWLn5m5k8n75CAR9MPcnR3SJWpdPc/O05uT/\n1a67enFLhnirsGEtMbBGSwz8+kXwN8HxN3Q49+xRWnx/++ixDQfrmdK6wUvBOKj+JtJ3OM9mwesP\n4lLZ4UfMWWedhd/vZ+rUqdx1113MmdNJoTzgnnvu4dJLL+Wkk07qtCXjVVddRV1dHVde2dLT6rrr\nruOmm26KLO6GKSgoYPHixVx00UVMmzaNyy+/PNope01yu3pc3qjRNWHCi2hr9tVx/jQtimPVnhqk\nhLljOu+tqWX0EnmIQXsoQYv1j8UAu4Vq5YPtnPr9gIS8UT0/x4gTtYzfVY/AsdfELnrmrII1T2oR\nNQOO7IOmzu0lN6PdjLJoMow7B758hInH3YDRIFi3v45TxxcCmn//cIOHm0flwuontMqiRR2n/5MG\nZ2GzGPmqpJZzp2qzkPIGD6V1TVw/t7hlx/wxEGiGxkOQMzzyQVTn8mJPS+rHu89JS0vj3Xff7bDd\n6WyJ/Lvkkku45JJLAFiwYAELFnTs9Xzddddx3XXXRX7+7LPPuOSSS8jJabELF198MRdf3NJoZ8WK\nFZHvzz77bM4+u8X1dzRIcsXv61TxTxiUSYbFyJpWquqz3dVkppmY2kXbvBybhXEDM9ss8K4/UI9B\nwLROWt9lp5vx+oNH1M0r5QiHJx6hIW6DEHDCD7WIl13vxd7vkz9ptfbn/fyIL1HrijGjPPl/wVOP\n/au/M2VINqtarQWt3KUV/jsr+BnUH4DZN0U9t9loYMaIXL5sdeza/dq9NnNEqxllbuhDoFZzDeSp\n7PCE4pZbbuHOO+/krrvuivdQ2tCdZutPCyEqhRBbWm27RwhxSAixIfR1ToxjzxJC7BRC7BZC3NmX\nA+8Kjy9Aky/Qplxue0xGA8eNyGVNaBFNSsnKXVXMHpXXqX8/zOziPNbtr8MXKqj15d4aJgzKIqMT\npZUdCi1tbPIdya+TWoRi4Xul+EHrYpU3Cj68BwJRMmgrtsLaZ2Dm9VqRtyOkzu3t6OoBGHKc5mpa\n9TBnDXaz8WB9JIN35a4qjskzk7/2z1p44fhzY57/lLEFfFPhjHTVWhNqrj6xddBBXsjw12mGPzcj\nnCui1pESgYceeojdu3czdmwP16uOEt1R/M8CZ0XZ/hcp5fTQ19L2bwohjMDDwNnAROBKIUS/pabV\nuzXDmtOJ4geYe0w+O8odHKx1s7G0gdK6Js4IlWHuihNGD6DJF2D13lrqXF7W7qvltNCUPhZhw9+g\nDH9savdq8fjdDKuMidEM838DVTtg1UNt3/N54M3vgS0P5v2iZ8OMpfgB5t8NJitXHvoDBH2s2lND\nQ5OPlbuq+XXWEs1Qz/9Np3X3z56iNXB/Z3MZwaDkg20VnDh6AObWoiRriNZ8PaT4cyLCInYRwURG\nLUp3TV/8jbo0/FLKT4GeLC3PAnZLKfdKKb3AP4GODrGjRF37ULsYnDtVe7heX1fKy6v3YzEZOHNy\n9wz/vHGFZKebeWn1ft5YX0pQwhmTOj82Sxn+rqndqy3s9kVY24TztEibj38PO5dp23xN8MZ3tQSh\n8x+CjAFHfNpAUOLw+CMf5B3ILILz/kp29df8Lf1J3lhbwpvrS5kf/IKTyp/XZgTHdJ6VOSQnnWOH\n5/DG+lI+31PN4QYP54Q+DCIYjJA7IqL4w/dXo0d/95fVaqWmpkYZ/06QUlJTU4PV2rH+2JHQm9Wf\nHwohrgHWAj+VUravKjUEONjq51Jgdi+ud0S0GP7Yrh6Aobk25k8YyKLluwlKybUnjIz9MLfDajZy\nzQkjeOjj3Xy4vYI5o/KY3MXaQLaOH8x+o3YPDJ3VN+cSAs77u+ZP/+eVMPo0qNkFdfvgrD/BuJ4t\nooXLcmdaO3mEJl8MtSX898e/Y+LunRzcVch1ls3IYbPhvx+MfVwrbjxpFN9/aT3feXYNRVnWjoYf\nND9/SPGHx6NHV+LQoUMpLS2lqqrzfhepjtVqZejQob06R08N/6PA7wAZen0Q+E67faLJtZgf5UKI\nG4EbAYYPj14q4UhoCLl6srsw/AC/WTCJ2pebyUo3c+vpR+aL+8Gpx/BNhYPD9R5+f8GULvdXrp4u\n8Hu1UM6pV/TdOdNz4Nq3tdLJuz/SGqz895+7VNydETas0cqBtOHkn+HNn4jvrT8ywtdI/fRbyDnz\nF90uNX3W5CKunzuSlbuqufu8iS05A63JK4aDq0FK0kxGrGYDjR79uXrMZnMkY1VxdOmR4ZdSRqoR\nCSGeAJZE2a0UGNbq56HA4Sj7hc+5GFgMMHPmzF7P9RyhGz9WBm1rhuSk86+b5/boOlazkcev7kZr\nwBBZIUUW/mBStKPxEMiglpHal6Rlwpl/0L76gJb7q+tHyDLxHMZOjBr/0CVCCO4+b1LnO+UWQ3Oj\n1hw+YwBZVrO6vxSd0qNwTiFE6/nmhcCWKLutAcYIIYqFEBbgCuCtnlyvJ4RdKZ1OxeNAiw9Wf4qs\nX2gIeQdzhnW+X5xxRO6v7rkFjyrtInuy083KlajolO6Ec/4DWAWME0KUCiG+C9wnhNgshNgEnArc\nGtp3sBBiKYCU0g/8EHgP2A68KqXsWMXoKBFWZImWxGI2GsiwGJWrJxb1IcOfndiGv/EIZpRHnZwR\n2mv9AUATF8rwKzqjS6sopbwyyuanYux7GDin1c9LgQ6hnv2Bw+PHZjF2Kx6/v8lKNyvDH4uw4s/u\n3eLV0caRSDPK8N+qoRTQ3E+q2Y+iMxLPKvYRDo8vMR7KKGSnm3UZddEv1B8EexGYOjbPSSTC/7+E\nuMesWVqFz7DhV4pf0QVJa/idzf7E8L9GQSn+Tqjfn/D+fWhxJSbMPZY9NDJbyrIqYaHonKQ1/A6P\nPzHUWBSyrGa1uBuLhoN9H9FzFHA0+7GaDVhMCfIItTb86SYaPX6VCKWISYLctX2P5upJEDXWDuXq\niUEwCA2HEn5hFzRXT0Is7IbJHtrKx28mEJSqNLMiJkls+BNX8WdaTZHFQUUrnOVaXXyduHoS6v7K\nGQZNddDsbAkZVuJCEYOkNfyNHn+3kmvigT3NhLNZTcU7EAnlTHxXT2OizSjDs6TGQ5GZiFrgVcQi\naQ1/Irt67FYTQQkeXzDeQ0ksdJK8BSFh0c2aTv1CJKTzIFnp4Xo9ah1JEZ2kNPxef5Bmf5DMBEve\nChNOKnM0K0XWhlACkh58/AkXLhw2/PUHI/0gwj0AFIr2JKXhT6jkmiiEDb9TRfa0pfGQVoc/zR7v\nkXSJI9FcifYiEEZoKI0IHocy/IoYJKnhT7AY63ZEDL96MNvSWKY1FtEBCRfVYzRB1mBoKFWKX9El\nSW74E0iRtcJuVYY/Ko7DkBml3nyCEXElJtr9FYrlD99fyvArYpGkhl9z9dgT7cEMoVw9MWgsg6zE\nN/wJVZmzNVmDofEwGZaQq0fdX4oYJKXhT6jKiVFQrp4oBPzgqoTMwfEeSZdE7q/0BBMWmYPAUY5R\nQLrZqBS/IiZJafi71RYvjqipeBScFVoDFh0ofmek5HeCCYvMQeBvAk89dqtJCQtFTJLS8CfsVDyE\nXUVddMRRpr3qQPG7vNr/LSMtShvEeJJZpL06yiNJggpFNJLU8Ce24k8zGTAZhPLxt6Yx1JVTB4rf\nHTb8lgS7v7JCH5qNh7GnmdSMUhGTpDT8zmY/aSYD5gRswgJaH1W7VT2YbdCT4m/Wip8lsuLPSDMq\nxa+ISWJaxl7iavYnXMvF9tjTTMrV05rGw2Awg21AvEfSJeEPbFuiKf5wKKzjMPY0M85mVZ1TEZ2k\nNPxubwBboqmxdtjTTMrV0xpHmWa4DIl/S4bLHWckmrgwp4M1J+TjN+JUJUEUMehOs/WnhRCVQogt\nrbbdL4TYIYTYJIR4UwiRE+PYfaGm7BuEEGv7cuCd4Wr2J57/tR1q8a0djYd14d8HcEcUfwKKi6zB\nIVePKeKSUija0x159SxwVrttHwCTpZRTgW+An3dy/KlSyulSypk9G+KR4/YGEvOhbIXy8bcjrPh1\ngMsbwJKoa0iZRdrirgrnVHRCl3eulPJToLbdtvellOG76ktg6FEYW49xef2JNw1vh/Lxt0LKUNZu\n4i/sghbVk5GowiJTU/x2iwmvP4jXr0p/KzrSF5LlO8C7Md6TwPtCiHVCiBv74FrdQjeuHuXj12hu\nBJ9LP4q/OZB4C7thMovAWUFmmgBUkqAiOr26e4UQvwT8wEsxdpkrpTwshCgEPhBC7AjNIKKd60bg\nRoDhw3vXgcnVrI/FXfVQhmgMhXLqRPG7mv2JF8oZJmsQyAB5NAJaaHNuhiXOg1IkGj1W/EKIa4Fz\ngatkjB6CUsrDoddK4E1gVqze06ioAAAgAElEQVTzSSkXSylnSilnFhQU9HRYQHgqnqCKLITdasLl\nDRAIqvaLOELJW3pR/InsSgz9DQcEqwFVD0oRnR4ZfiHEWcAdwPlSSneMfTKEEJnh74EzgC3R9u1r\nXDoJ54SW9P+UxlGuvYYTkBIctzeQuMIiZPiz/TWAMvyK6HQnnPMfwCpgnBCiVAjxXWARkInmvtkg\nhHgstO9gIcTS0KEDgc+EEBuBr4B3pJTLjspv0QpfQFvQStgHM4QqzdwKZ4X2qhPD72r2J27UWMjw\nZ/qqAGX4FdHp0jpKKa+MsvmpGPseBs4Jfb8XmNar0fUAdyi5JmEfzBC2kOEPjzelcVaCJRMsGfEe\nSbdwewOJ6+qxF4IwYPdWA6PVOpIiKgkYiNw7wjd6opdsCIcDupWrR1P89sJ4j6LbuL0JrPgNRrAP\nxNqkzaLUjFIRjaQz/GFDaktwwx8OB1TZlWiK3z4w3qPoNs7mBF7cBbAPxOJRi7uK2CSd4Y9UTkxU\nRRYiHA6oFD+6UvyBoMTjS/A1JPtATG7Nx6+EhSIayWf4vQlaObEdEcWvfPzgqNCN4ncnahOW1tgL\nEa5K0kwGJSwUUUk6w+9O1Frp7Ygo/lSfivuaoLlBN4q/JXgggYWFfSA4K7FbDMrVo4hK0hl+pfh1\nhrNSe9WJ4g8HDyS0sLAPBBlgcJpbRY0popJ8hj+k+BM9qiccFdKU6lNxnRl+fSh+bfY0xORQil8R\nlaQz/C1RPQmsyACz0YDFaFCKP5y8pRNXT9iQJnTwQOhDdJCxUcXxK6KSdIY/rPht5gR+MEPY0ozK\nxx8x/HpR/GFXT+Ir/oHGBiUsFFFJOsPv9mqN1k2J2CSjHRkWk3ownZWAgIz8eI+kWyRso/XWhD5E\nC6lXil8RlcS3jkdIQldObIfNYlThds4KrcG60RzvkXQLtx6CB9LsYM5gAPVqRqmIStIZfndz4rdd\nDGNTfVF1l7XbkiCYwIYfwF5IbrBOLe4qopJ0ht/Z7E/4iJ4wGUrx6yprF1oUf3qiiwv7QLIDdbi9\nAWK0y1CkMEln+PXQaD2MzaIUv94Uv7M5gMVowGJK8EfHXkimvxZ/UNKs+u4q2pHgd++Roycff0Za\niit+KXWp+BN6YTeMfSA2n9aMRSVxKdqTdIZfVz7+VI/q8TRAoFlXij+hG623xj4Qq68BCz4V2aPo\nQNIZfpcO+u2GybCkeBx/OGtXJ523QE+KX5tFDaBRtfdUdCDpDL9bB/12w9jSTLh9Kbz4prOsXQj1\nc9aDsAjNogqEiuVXdCTpDH/CN8lohc1iRErw+FJ08U1nWbugVVPVk+IvEPU4Uz2AQNGBbhl+IcTT\nQohKIcSWVtvyhBAfCCF2hV5zYxx7bWifXUKIa/tq4NHQS6P1MOF6Lyk7FY8UaFOKv8+JKP6G1HYn\nKqLSXcX/LHBWu213Ah9JKccAH4V+boMQIg+4G5gNzALujvUB0RfopdF6mLABcaeqInNWgNEC1px4\nj6TbuPSSJ5JRAEAB9SqJS9GBbhl+KeWnQG27zQuA50LfPwdcEOXQM4EPpJS1Uso64AM6foD0Gboo\noNWKsMsgpRW/fSAIEe+RdJuEbrTeGpOFYHqepvhTOXJMEZXe+PgHSinLAEKv0ebrQ4CDrX4uDW3r\ngBDiRiHEWiHE2qqqqh4NKFKZUw8PJq0Uf8oafn3F8IN2j+lFWGAvDPn4U/T+UsTkaC/uRpNyUUNY\npJSLpZQzpZQzCwoKenSxiOLXgw+WVoo/ZV09+sraDQQlTT795IkI+0AKRUPqCgtFTHpj+CuEEIMA\nQq+VUfYpBYa1+nkocLgX1+yUSJMMnSgypfj1pfibfDop0BZC2AdSaGhIXWGhiElvDP9bQDhK51rg\nP1H2eQ84QwiRG1rUPSO07aigl0brYcIGJCUfzGAA3NW6Uvzh6Bi95IlgLySfelweX7xHokgwuhvO\n+Q9gFTBOCFEqhPgu8CfgdCHELuD00M8IIWYKIZ4EkFLWAr8D1oS+fhvadlTQS6P1MGEDkpKK31UF\nMqgrxR+eUeoiqgfAPpB0mvE3O+I9EkWC0a07WEp5ZYy3Touy71pgYaufnwae7tHojpBw9ILuFH8q\nRl3oMXlLD43WWxP625rcPQuWUCQvSZW5G05N18uDaTUbEILUTLCJJG/px/C79NBovTWh2ZTVUx3n\ngSgSjaQy/BHFr5MHUwiRun13dVinJ6L4deTqAbA2K8WvaEtSGX5Xs34arYdJtxhTM8EmbPgz9GP4\nXV69KX7N8Gf4jtqymkKn6MdCdgM9NWEJk7LtF52VkJYFFlu8R9JtwlFjulH86bkEMJLpV4Zf0Zak\nMvx6asISJmXbL+oshh9aRfXoZA0JgwG3OY/sgDL8irYkleHXUxOWMCnbftFZqSs3D+io0Xor3JYB\n5Ml6/IEULf2tiEpSGX49NWEJk7LtF50VkKmfiB7Qwm510Wi9Fc3WfApEPW5fCt5jipjo5w7uBq5m\nnSr+lAznrNJVKCdoYbd6Exa+9ALyRaPqwqVoQ5IZ/oBukrfC2Cym1Ivq8TVBc0OkZrxecHkDuhMW\nAVsB+TTg8njjPRRFApFchl+PPn6LMfXq8esweQt0VIu/FTKjELMI4GmsifdQFAlEUhl+Xfr400yp\n14FLp4bfqada/GEyiwDwN5TFeSCKrvj6QB3LtpT3y7WSyvDr0sdvMeIN9QpOGVz667ULOmq03gpT\nlvbhGnBUxHkkiq54de1B7vrPlq537AP0ZSW74L2fnKw7RRauK9TkDegqWqRX6LBcA2g+/twMS7yH\ncURYsjXDL5XhT3hczYF+ywrXl5XsgpH5GfEewhETVpBun59szHEeTT8RdvXobHHX7fXrp1xDiLSc\nwQAYXNH6JCkSCVdz/1UeSBGJmbikp2IzFmcl2AaAUV8fdK7mgH7KNYSwZebQJC0YVWnmhKc/g1OU\n4Y8zYQWZUtm7zgrdLeyCPhW/zWKiSmZj8SjDn+j0Z3CKMvxxJuzjd6ZSgo2zUndunmBQag+mzoIH\nDAZBjchVNfl1QH8GpyjDH2fCbfxSKqRTh4o/XPJAN20XW1FnyCVdlWZOeDRhoRR/ShCe2qVMEpeU\nWr9dnUX06K7ReisajbnYfSqBK9HRxeKuEGKcEGJDq69GIcRP2u0zTwjR0GqfX/d+yMlFeGqXMmUb\nvE7wuXWn+F2R7m76U/xO8wDsgQYI+OI9FEUMpJS4+lHx9/gullLuBKYDCCGMwCHgzSi7rpRSntvT\n6yQ74XDOlCmi5dRn8lZLP2f9KX6XeQA0oc20sgbHeziKKDT7gwSCMvEVfztOA/ZIKff30flSBluq\nhXPqNXkr3Ghdhz5+T9oA7RunSuJKVPq7X3hfGf4rgH/EeO8EIcRGIcS7QohJsU4ghLhRCLFWCLG2\nqip1Qs+MBoHVbEidcE6d1umJPJg6NPxea772jVMlcSUqkRmlXhS/EMICnA+8FuXt9cAIKeU04CHg\n37HOI6VcLKWcKaWcWVCgr1C/3mJPM6VOOGcka1dnil9vjdZb4U8P/a2V4k9Y3P28htQXiv9sYL2U\nssNdJaVslFI6Q98vBcxCiPw+uGZSkVI1+Z0VIIxgy4v3SI6I/lZkfUkgI6z4leFPVMLCQk8JXFcS\nw80jhCgSQojQ97NC11NxZe2wWYwppPgrtOQtg76Uc3gNRjeN1luRnp5Bg7SpQm0JTDiPp78Uf6+u\nIoSwAacD32u17SYAKeVjwCXA94UQfrS4giuklLI310xG7Gmm1PHx6zCGH1pKaugxjl8r25BDhqMi\nuaoyJhHOSPBAgodzAkgp3cCAdtsea/X9ImBRb66RCtjSTDQ2pUiMtbNCl4bf2ayVzTYb9ZfzaE8z\nUiVzGO5Qi7uJijuyhqQfH7+il2RYjKkVx6+ziB7QZ4G2MBlpJqrIRqionoQlnCCoJx+/opdkpKXI\n4q6UIcOvR8Xv112BtjBhV4/BrQx/ohIuCaIUfwqRMg3Xm+og6NOn4m8O6LJAG2hrSFUyG6PPCV5X\nvIejiEJY8aebleJPGWxpptRw9ei08xZo4XZ6XNgFzX1QTbb2g3L3JCTuZj82ixGDQfTL9ZThTwDs\naSZ8AZn8Ddcj5Rr0p/hdzX6dK/4c7Qdl+BMSVz/3elCGPwGwpUoXLleoFIcODX9/1krva2wWI1Uy\nrPhVLH8iogmL/ru/lOFPADJSpQuXTgu0gfa/0WNJZmiv+JXhT0Tc3v4NHlCGPwEIF/5K+sgeZwUY\nLWDNjvdIjhi3N6DLAm2gRfXUkEUQg3L1JCiu5kC/JW+BMvwJgS1VavKHY/hF/yxg9SWuZv0u7lpM\nBkxGE25zDriU4U9ElOJPQTJSpSa/TmP4/YEgzf6gbl09oJUCcBjzlOJPUFxepfhTjoxU6burwybr\n0Krtok5dPaC5exqMecrHn6C4+zlBUBn+BKCl726SG35HGWQWxXsUR0yk+5ZOo3pAW+CtNeQoxZ+g\naMEDSvGnFGHfsTOZXT3+ZnDXQOageI/kiIkU0NKz4k8zUkuOpvhVgdyEQkrZ78EDyvAnAOHEIHcy\nL+6GXQy6VPxhV4++FX+lzIaAFzz18R6OohXeQBB/PzZaB2X4EwKryYgQLb7kpMRRrr3qUPFHum/p\neHHXZjFSEVBlGxKRcBOW/kwQVIY/ATAYBDZzkpdmdpRpr3pU/KEPZL2WbADNTXXYn6X9oBZ4EwpX\nP9fiB2X4E4aMZO/CFW77Z9eh4Y8ofn27ekr9mdoPqgVjQuHu51r8oAx/wpCRZkruOH5HGRhMYBvQ\n9b4JRliR6Vnx2ywm9nlDrh7H4fgORtEGVz/X4oc+MPxCiH1CiM1CiA1CiLVR3hdCiL8LIXYLITYJ\nIY7r7TWTEVuyd+FylGtq36A/rRHxwerY8NvTjNQFrEiLHRrL4j0cRStccfDx99WdfKqUsjrGe2cD\nY0Jfs4FHQ6+KVmSkmZI7gUunMfzQUjzP1k9NMo4G4YXpoL0Io0MZ/kTCFYdw4f6QXwuA56XGl0CO\nEEJ/oR1HmQyLMbmLtDnKdWv4tToq/dck42gQdlP5bANbFtoVCUE88kT6wvBL4H0hxDohxI1R3h8C\nHGz1c2lom6IVtjRTcpdl1rXi798mGUeD8MKh1zZQuXoSDL26euZKKQ8LIQqBD4QQO6SUn7Z6P5pM\n6pA6GPrQuBFg+PDhfTAsfWG3mCK+5KTD59GShnRq+N3e/m2ScTSIlP5OKyTLUQbBoC7XW5KRsODr\nz+CBXv/npZSHQ6+VwJvArHa7lALDWv08FOgQViClXCylnCmlnFlQoL+erL3FlpbEDded+k3eAk2R\n6V3xh42KK61Qa3jvronziBRhnB4/BqGjBC4hRIYQIjP8PXAGsKXdbm8B14Sie+YADVJKNddsR4ZF\na7guk7GOSiRrV5+K39Xs13W5BmgxKo2mfG2D8vMnDM5QP2fRj30qeitjBgJvhgZsAl6WUi4TQtwE\nIKV8DFgKnAPsBtzA9b28ZlKSkWYiKKHZH8Sq4+iRqESydvWp+J3NfvLtlngPo1eEFX99a8M/aGoc\nR6QI4/D4ybSa+/WavTL8Usq9wLQo2x9r9b0EftCb66QCYR+yw+NPQsOvb1ePs9lPcX5GvIfRK8Ku\nqlpDKIGuUSVxJQrOZl+/Jweq1Z0EwW5N4obrjjKt1256brxH0iMcHl/k/6NXwoalmhxAKFdPAqEp\nfmX4U5LMNG2q5/D44jySo0BjKJRTh712ARrj8GD2NVazAYMAp19o7S+V4k8YnM3+fhcWyvAnCGHD\n4vAkoeJvPARZQ+M9ih7R7A/g9QfJ1HG5BgAhBBmWUK5IZlGL+00Rd5wev3L1pCrhxZ2kVPwNpZCt\nz5w9Z+iDuL8X344GtjSjliuSOVi5ehIIR7Ny9aQsSav4g0HNrZClT8PviBh+fSt+0CLHnF4/ZA1S\nrp4EwuHx9buwUIY/QUhaw++q0hKGsvXp6olHVuXRwp5m0tp7Zg6Gploto1oRV3yBIB5fULl6UpXw\nPz7pDH9jqfaaNTi+4+ghjSHXW1K4eixGrS5MOJFOuXvijitOwkIZ/gTBZDRgsxiTz8ffcEh71amr\nx5lErh57uPR3ePbVeCi+A1JEhJ6K6klhMq2mJFT8IeOiU1dPMvn4baGyIGSHSmfVH+z8AMVRJ+xK\nzFKGP3Wxp5lwNCeZ4m88BCarLlsuQkuUVTK4ejLSTDibAy0RVg2l8R2QokXxp6nF3ZQl02pOPsXf\ncEjz7+s0eSuZFne1Zj9+MKdDRgE0HIj3kFIeZ0joKVdPCpO0rh6d+vdBU2RpJgMWk/4flYw0E25v\ngGBQau4epfjjToviV4Y/ZcmympNzcVen/n0Il2vQv5sHWoyL2xfQ/ifKxx93lI9fkXyKPxjQQgZ1\nrPidcciqPFqE2y+6mv2QM1xT/MnY/0FHqKgehba4m0yG31EOMqDbcg0QzqpMDsMf6cLVHArp9Dep\nTlxxxunxYzQI0vu5FLsy/AlEptVMky+APxCM91D6hvr92mvOiPiOoxfEo4DW0SJck9/VHGgJ6WxQ\n7p54Eo/uW6AMf0KRmWw1+ev2aa+5I+M5il4Rj1rpR4tw+8g2SVzKzx9XGj3934QFlOFPKJKuXk/d\nfkDoenE3HgW0jhYZllaunpzh2kYV2RNXGpv8ZKX3//2lDH8CETYwjckS2VO3T1vYNaXFeyQ9xtGc\nPK6esIFp9Pi0bmhmmzL8caahyUuOngy/EGKYEGK5EGK7EGKrEOLHUfaZJ4RoEEJsCH39unfDTW6S\nTvHX74dc/fr3A0GJs9nf76F2R4uwgal3+7SEuuxhLeswirjQ0OQjOw6Gvzd3tB/4qZRyvRAiE1gn\nhPhASrmt3X4rpZTn9uI6KUNWWPE3JZHiH/1f8R5Fj3F4fEgJOTZLvIfSJ4QVf0P4/sorblmHUcSF\nerePHJuOFL+UskxKuT70vQPYDug3bi8BCN8A9clg+H0eLYZfxxE9dW7t/5CbkRw+fqNBkGU1aYof\nIG8U1O5VsfxxpD5Oir9PfPxCiJHAscDqKG+fIITYKIR4VwgxqS+ul6xEDL/bG+eR9AH1oTowOnb1\nhP8POenJofgBsm3mVop/FPjc4KyI76BSFI9P6+ecrSfFH0YIYQfeAH4ipWxs9/Z6YISUchrwEPDv\nTs5zoxBirRBibVVVVW+HpUvsaSZMBtGiyPRM2Hes41DO8MwrHg/m0SIn3dIiLPKKtdfavfEbUAoT\nfs7jISx6ZfiFEGY0o/+SlPJf7d+XUjZKKZ2h75cCZiFEfrRzSSkXSylnSilnFhQU9GZYukUIQY7N\nHHEx6Jqw71jHrp6GyIOZRIbfZm5xJeaN0l6V4Y8L4ZmXrlw9Qks1ewrYLqX8c4x9ikL7IYSYFbqe\nyhHvhBybhYamJHD11O3T6vDbB8Z7JD2mLuzqSZLFXdCMTPgDjezhYDApwx8nIq7EOMwoexPVMxe4\nGtgshNgQ2vYLYDiAlPIx4BLg+0IIP9AEXCGlWknqjJx0M3WuJFD8tSWa2jfoN1UkPBWPhyI7WrRR\n/EaTlsilDH9cqI+j4u+x4ZdSfgZ0WmBCSrkIWNTTa6QiOTYLpXXueA+j91R/AwXj4j2KXtHQ5CPL\nasJo0GcTmWiEffzBoMRgEC2RPYp+R5euHsXRIad11IVeCfigrgTyx8Z7JL2i3u1NKjcPaPdXUILT\nG0oSzBulzc7URLzfCbvcdBnVo+hbcm3miG9Zt9Ttg6Bf94a/Lk7JNUeTSBJX61j+5kZVnjkONDT5\nMBoEmapImyLHZsHjC+LxBeI9lJ5T/Y32qnPDX9/kSz7F37psA0DeaO21ZnecRpS61Li0Oj39XZIZ\nlOFPOFqSuLp29wSDkv01Lnw9qd/va4KyjVC1U+uU1ZdEDP8xfXveXiKl5GCtu9sfqg3u+BTQOprk\nZWgfZDWuZm1D4XjttXJ7h32llJRUu2jy9uD+CAahYhuUbQJ/c0+HqzsO1rqpdXVvxl7jbCbfHp8C\nhslRfSqJyA0pzDq3l6Jsa8z9PL4A1z+zhlV7azim0M7LC2dTmBV7/wg+D3x6P3y1WJviA+QWw5l/\nhPHn9MWvANW7tDBOa3bfnK8P8AeC/PDlr1m2tZwhOem8fMNsRgzI6PSYGqc3YiiThYJMzdBUO0PG\nKXsYWOwdDL+Ukl/+ewsvrz5Ars3M89+ZzZSh3fx/bvkXfPDrliYvadlwwg/gpNvAmFwfpK158P2d\nPPTxbjIsRp64diYnjo6ashSh2tlMfmZ87i+l+BOMAWFF5uxcNTy6Yg+r9tZw/dyRHKpr4qevbez6\n5I5yeHI+rHwAjpkPlz4H5y8CSwb880r4oo8CsKp3JZyb55EVe1i2tZyr54yg0ePjjjc2dbq/xxfA\n0eyPGMpkIawwqxwhFS4EFIyHyra1Fd/bWs7Lqw9w8XFDsVlMfP+ldbi93agau+JP8Pr1kJEPFzwK\nlzwDxSfBij/C02eCMzmz8pfvqOShj3dz3rTBFGVb+emrG7v8e1U7vXFT/MrwJxhhQ1Pp8MTcx+ML\n8PTnJZw1qYi7z5vEz84cx8pd1XzyTScPVf0B7cGr3QvffhUufQYmXQDHXQ0LP4KJC+D9X8Lm13v3\nCwSDmnpMoFDOOpeXxz/Zw1mTivjdBZO57fSxfLm3lnX762IeEzaMyWb4M9JM2CxGqp2t3C+FE6Bq\nR+RHKSUPL9/DqIIM/u/iKfz5smmU1jXx7Bf7Oj/51y/Binth+lXw3Q9h+rdh8kVwxUuayKjYBs+c\nrQmQJGPR8t0MzU3nz5dN4/cXTKGswcO/vz7c6THVytWTAASDWmPw/piKBoOwdzlsf0vzgXrqIT0P\nhs2iaPxFQCtFFoUPtlXg8Pi5+gStHMLVc0bwzOclLPp4F6eMjVLuoqkeXrwEmurg2rfxDzqWpz7Z\nwzuby8i3p/Gj08Yw/aInobEMltwKQ2d2q8aOlLLjwlT9PvA6oGjqkR/bBZUOD/cv28mWw40cPzKX\nn54+rluhcM98XoLLG+DW07VZyGUzh/HAezt5fV0pM0bkRj2mypmchh801d/B8H/9gqbG7QVsL3Ow\n+VADv10wCZPRwOxRAzh5bAFPf1bC9ScWk26J0hi8cju8cxsUnwznP8R/NpXzwqr9GA2C73yrmDMn\nXaC5/168GF66BK5/F9IyOx3nEd0jUkLpGtj0Khxer93z6TkwaBpMvEAbVxfncnv9/O3DXazcVc3o\nQjs/O2Nsl+5AgHX7NRFxz3kTMQvJnBGZjC/K5J9rDvDt2cNjXsvtDTDArlw9/YvXBRtehn9eBQ+M\nhd/mwu/y4YFx8I9vw9cvgvcoJFJVbINnz4EXL9LUtTUbBh+rlTdY+zS2Z+fzfNp9+KtiR1n8++tD\nDMq2MmfUAAAsJgPXnTiSNfvq2HKooe3Ofi+88j+a0r/8JYKDj+P2NzZx77s7MBkEmw81cNnjq/jg\nmzq4+EntmLd/EjOu2+Hx8ev/bGHKPe8x5Z73+eWbm9tWEy3frL0WTelwrD8QZPGne5j9xw8Z96tl\nLHxuDXurnN36sx2ocXP+Q5/z1sbD5GWYeXn1AS5fvKrLnAe318+zX+zjrElFjCvSDE1GmolTxxfy\nwbZyAsHov2dE8cdJkR1N8u2WtsKicIL2WqX5+d/bWo4QcM6UQZFdbp43mmqnlzfWR+nYJSUsuU3r\n6HXxUzz5+X5+/M8NNHp8VDmb+d4L63j8kz0w4gS47HntGXj1Gi3fox3+QJBFH+/i+D98yLi7lnHT\nC+vYX+Pq/Beq2AbPL4CnTse//kW+LvfydlUhGysDBDa8As+fD4+fDLs+iHmKRo+PKxd/yeKVe8mx\nmVmxo5ILH/mCPV3dn81Otry9iOesD3Lt5/Pht3mI3xfyuvu73FTxG+pXPacFUrSj2qE9M8rV0194\nGmH5HzUD/+/vQ/kmrVnIKXfAKXdq35dvgv/8AP4ySfNZNnfPOHWK1wUf3A2Pn6RF0pz3d/jfPXDN\nv+GSp+H6d+Bn38Dpv2WG2MmNW6+Cdc92MMC+QJBVe2s4feLANhmllx0/DJvFyNOfl7TsLCW8/WPY\ntxIWPAzFJ/Hi6v38a/0hbp0/ln/dPJf3f3IyE4oy+dE/vmZbUy7816+02cjOpR1+hRpnMxc98gUv\nrT7A/AkDOXNSEa+sOciFj3zBwdrQh2T5ZhBGKJzY5lh/IMgt//iaPy7dwdiBmVw1ZzirS2pZ8PDn\nrN1X2+mfzu31c+MLa2nyBXjz5rm8tHAOT193PHuqnNzZha/+X+sP0ejxc8PJxW22nzW5iGqnN+a1\nk9XVA9rv1Fbxh/5XFVsBeH9bBTNH5LYxSrOL85g6NJunPysh2P7DcvNrcOALmH83u1zp3LdsJ6dP\nHMi7Pz6Z935yMudOHcS97+5g6eYyGDMfzvsb7PkYlrQVGL5AkO+9sI4H3v+GqUOyueL4YXy+u5rz\nF33Ouv0x7pGvX4LF86B8EytH/5Sprof5Vfa9bDnxL9xsvJvJrodZO/330OzQZhqvXQeOtmWog0HJ\nT1/dyNbDjSy+eiYv3zCHt2/5FgDff3Fd9KgmTyN8/HuCD4zj2uoHmWatQIw9S7Mhp9xOcOS3ONaw\nm5z3fgR/ngjL79WOCRGZUcZLWEgpE+5rxowZss/xNkn5+UNS/mmklHdnSfnK1VLu+1zKYLDjvsGg\n9t7LV2j73j9Wyq9fkjIQ6Nm1dy6T8s+TtXO9ebOUzupOd1+46G258Y/ztP3/dZOUXnfkvbX7auSI\nO5bIpZsOdzjurn9vlmN+sVRWNDZpG5bfq51j+Z+klFIerHXJiXe9K69+arUMtvq9Kxqa5Ow/fChP\nvPcjWdvokvLhOVL+bbqUfl9kH1ezT56/6DM59pdL5We7qtqMZ+o978kT7/1IHqhxSfnSZdrxrfAH\ngvKWl9fLEXcskU98uidy7YO1Lnnq/cvllLuXyV0VjVH/FsGgduzIO5fI5Tsq2rz38PJdcsQdS+SS\njR3/FuFjT3twhTzvoRTwTQgAABjVSURBVJVtfl8ppXR4fHLUz9+R9y3bHvXYP7+/U468c4n0+nv4\nP09gfvGvTXL6b95ru/GBcVK+vlCW1TfJEXcskY+t2N3huH9/XSpH3LFEfritvGWjt0k79vF50u/3\nywse/kxO/817srLRE9nF4/PLCx/+TE646125szz0f/7o99q9+cn9Ukrtf3XbKxvkiDuWyOe+KIkc\nu7/aJefdv1xOvnuZ3HKovuW6fq+US2/XzvHsufLZ99fIEXcskT/+x3rpD2j/a6fHJ6964ks58s4l\n8q31JVKuuE/K3+ZLee8wKde/GHn2w/fRUyv3tvl9V+yslCPuWCIfeG9Hq+v6pFy9WMr/GyXl3Vly\n458XyMt/+RdZFX7mQgSDQTnvvo/lvY883mJH/m+UlF8+JqWvWb63pUyOuGOJ3HSwXvYVwFrZTRub\n/Io/4If1z8NDx2mLl4Onw40rtCnniBOj+/2E0N678h/wnfche4g2O3jiVNj/Rfev3XBIc7O8fBmY\n0+G6pXDBw5AxoNPDzDmDuM38K20WsvFlePqsSFPsVXu0DMvZozqe4/q5xfiCQZ7/Yr/mxgovtJ1y\nuxae9+YWAP544eQ2vtPCLCuPXz2DKkczP3p1M4FTf6W5hja9AmhK7Icvf83m0noWffs45h7TEqY2\nY0QeLy2cjbPZzxWLv8R/aGMbN08gKLn99U28tfEwd5w1noUnjYpce2iujee+MwuLyci1T6+hsrHj\ngvaTK0t4a+NhfnbGOOaNK2zz3o0njWLq0GzufmtLSyZqKz7bXc3uSifXnTiyg6/YnmZi6tDsyN+z\nPVXOZvJsFszG5HtECjLTqHP72uZ/DJkBh9axukT7e0QLRTxnyiAGZVt56rNWs8r1z2ud1ubfwzNf\n7OfrA/Xcfd6kNjOlNJORR/9nBjaLiZteXIfD44NTfwFTL4ePfwebXuP/lu3kjfWl3Dp/LNecMDJy\n7PABNl5cOBt7molrn/6K3ZUOcFXD8xfA6sdgzg94auSD3P1RBedNG8wDl06LzIQz0kw8cc1Mjh+Z\nx62vbePjgdfA97/QZjj/uRlevJjVX2/kgfd2cu7UQVw/t+W6AKeMLeDCY4fw+Cd7NZfP3k+0GfvS\nn0HBeOqueo9Lar7H6BmnkZ/ZNpRaCMGJx+TzQtlw/Je9BDd8rLnU3r0dFs0kc+frGAiqcM4+x10L\nXz4Gj8yGt26BzCK45i24+k3Np95dhs/WIhQuXAyuKi0q4ZWroXxL7GM8jfDxH+ChGZpf8bRfw02f\nwci53bpkYaaVCqdfeziu+AfU7IHHT4F9n7Nqbw3jizKjxpcX52dwxsSBHFz1OvKtW6D4FDj3ryAE\nSzeX88k3VfzszHEMzbV1OHbasBx+s2ASK3dVc19JMXLQdPj0PgI+Lz99dSMf76jkdxdM5vSJHcss\nTx6SzUsLZ2Px1GBylVGXpSUFNXkD3PrKBt5YX8ptp4/l+/NGdzh2WJ6NZ647njq3l/95anWb5Jfl\nOyq5993tnD25iJujHGsyGrj3oinUuX38aVnHBKQnVpaQb7fw31MHdXgP4IRRA9hU2oCruWPYXZWj\nOSndPNDiV24TMjxkBtTuYeM3+8hMMzFxcFaH48xGA9eeOJIv9tSw9XCDlpj12V9g+Ansz5rBA+/v\n5L/GF7Jg+uAOxw7MsrLo28eyv8bNba9uxBeUcP5DMOJbBN78PptW/oerZg/nR6d1TPobkpPOiwtn\nA4K7H38Z76Mnw6G1BC94nL8Yr+N37+7inClF/OWyaZjafVCnW4w8de1MJg7O4vsvrmdVQ54mwM6+\nn+D+VUz+zxn8KnsZ/3feqKgLyT8/ZzzHmfdS99Ql2lpBs1MTjdct4bFd2fgDQRaeNCrq33n2qAG4\nvAG2Hm7U/r7Xvg1Xaet6J2z6Je9Z7qBgx4vaQnQ/kzxRPVLCsp+Dq1IzlOWbQAa1P/jlL8L4c7tc\n1Y+JwQDTLocJ58GqRdrNvv0tGDYHxp4BA6eA2aop/H0rYeu/weeCSRfB/HuOuP1gQWYaDo8fjy+A\ndfw5mlr455XI589nrO9qOH5hzGPvGLKFYbsfpCJjPEWXvwAmC40eH795eyuTh2S1UVPtuXLWcDaV\nNvD4pyUUjr2C75bdyZOL/shbFcdzx1njuWp27N9j8pBsnp0fhA/h5k9N5FasY+PBBg7VN3HHWeOj\nGv0wU4Zm8+S1M7n+mTUsePgzbp0/lrIGD3/98BsmDMrigUunxYzumDQ4m4XfKubxT/dy4bFDmVWc\nB8DqvTV8+k0Vd549njRTlCgU4ITRA3hkxR7W7KvtMJuoTGLDXxRK9Ctv9LQkCQ6ZAYBz72qOLz4p\nZkXSK48fzt8/2sUjy/fw8Jh14DhMYMEj3P7GZswGA39oN5tszZxRA7jrvydwz9vbuO6Zr7h6zghW\nWu/k6sD3eT7tfsS4CQjRMSgAYHR+BktP3EHOp3dT48jisUEPsnrFYHaU7+Li44byfxdP6WD0w2Ra\nzTx3/Swue3wVC59bw0/PGIfVfBYv+G38QjzFdzzPw8P/gfHnaTkHmYPA74HyzRTufJd/spbGJhub\nxv6AqZfdBeZ0Kh0enlu1jwXTh1CcHz3yZ074XiypYdqwHM3+jDkdRp/GC888xOyDT2N692fw3s+1\nWXL+GE2gnv7bqOfrS5LH8AsBO5Zo4ZhZQ+Ckn8GEc7Vwrr7CYoNTbtcM74aXNHfKR+3+SWnZMPlC\nbZ8jmVm0ojBkcMobPIzMz4CCsXDDx9S/cC13H3qGyrJ9UHWftj2Mqxo+/j2j1j3DLttULq/7Ca82\nGim2aAtXNS4vT1wzs8sSw3+4YDJ5GWbu/QRmmEZxTv1L5F14DZfOjq5qWjPCtYmg0crQcXNYe8jB\nqIIM7r90apcZjKC5Fl6+YQ63vbqB217VktHmTyjkwcumk9FFEasfzx/DO5vL+N/XN/L6TSeSkWbk\n1//ZSlGWletOHBnzuJkj8jAbBav21nQw/Ifqmhg/vjDGkfpmaF46AKV1bqYPy9E2Dj4WiWCgYyvH\nnLgg5rHZNjM3njyKRz7chqf0AazD5vDXPYNYXbKH+y6ZyqDs9E6vfd1cLRz0t29v4/PdNZgMgkEn\nPs24sl8gXr0a5tysPWPprcJsK7bBB3dRuPtDfKNO4+XMn/J5iZ9Mq5G/Xj6dBdMHdxn2mZth4cWF\ns/nhy+v57RItWW3KkJGM+PY70LRNy2Tf/hZs+P/2zj26qurO45/fzesSQp6EEMiDhDeI8gghUMUH\nVim4wCooVkcqTi3o6JqO49QuXV2zdLkqth2dWeJYbatUl4LQUvE1FAWLlIBAxfAOSWhIIBDCI+EV\nCMmeP/a54RDuTe7NfZHc/VnrrnvOfpzzvb9zzu/us5/vXp4xczQt332en+wYyca9TSyva2JohpNn\nV+ygpQWemDLY4zn7JDrJ792TTRXHeWSyreDjcPB/agJ/Sh/Hiu8nwK4/Q9VmOFAMMT2N4/eZn7RT\n/RJI4lNh0uP6c+aYrg+/eE73U04bBA73JUxvyUnVVTEHjp/Vjh/AmcTiAS/SWPkKPz32ESwaD/3G\n6hGy547D/q+g+TxM/BcSxv8Hjte+ZvbrG8hM6sGumgZ+fscIXeroAIdDeOr2Ycz7Th713zaSvXoe\n2bHFQMeOnwPFOLLG8cv7Cjv1u8flprDmyZvYXdNAojOGnLQrq6TcER8bzX/PGc39v93EjFfX0yMm\nin8cO8NbDxXijPF8LXrERnFdVjKbKi7vMXLuQjN1p8+Tndq+E+uq9E92OX5bN0NnIqcSBzH+xF6S\n8tpvg1pw00Ci//4WzrOHef704/xubTmzxmUxe1yWV+e/d3wOd1zbj9Ijp8hN66mrLS+s1NM8bHwN\nti7WJW9nkh4fcLgE4hJh6ovEFP6YJx0OnuzE785IdPLBjydSfvQ0zS0wJCPB+sMo0GNXmpt0bcGZ\noxAdp5/l+FQcwHMjz3H3axu45/Vi8tJ7UlJdz7PTh3ss7buYkJ/KxyU1NLeoywpdB0+cY2T/JN3m\n2G90J36Nf3TfOv5Q0TMNssfrASLpQ/12+kCrwztw/PJxBMUVJ1ifcT/yxDbdbhAdB5V/0w2/Y+6H\nRzfB7S+QmZrE+z8q4rrsZBwOeGnWtcy7Ps/dqTySlhBH/qS79ECsr36lG8nbo7FBT/qWU+TTedoS\n5RCu6Z/ktdN34Wpkzk2LJ7FHDG8+WOB+MFsbivLT2H6wntO2en7XQjjZqb5p6Cr0csaQ1COGgycu\n719eEjuGQsceRvZpv8ExTlpYEP0R/3AOZ33zKP711sEsvPtanwbk9YyLZkxOyqW2qth4mP4r3fh6\n7Ww9tXflBu38v/s8PLENihb4vaKbiDCoTy+G9u11pd6oGD1pXd4NkF2oC3gW/ZN7sGz+RIoGptHU\nrHhu5kge9uKZmpCXxqnGi+yuudSVs6VFcai+kazk8BUsuleJv5uQ0ctJbLTjUt949DQN3xw4ydxJ\nuZCQDjc8qT8eGJzRi7cf6lzJuxUR3bNo6f2wYzlcN8dz2oq1eg7+gVP8O6cfjMtNZckjE33KMyE/\nlVfXlrG18kTrH0WV5fjdNYJ3F7JSerT+ThefnB7K9dIE1Ztg4M2eM5csJaqhigE/+IBVQ24MrLCM\nkbqf/1VIdmo8bz5Y4FMe1yDL4vJjXNNfT3JXe+o8Fy62kJUSPsdvSvxXIQ6HkJ3Sg8pjlx7MrZUn\nuNDcwsSB7b+GB5yh0yDjGj2jZ3vTN5f+RbdvZE8InbYAMC43hWiHsLHiUrdOl91zummJH2BA755U\nHL00IvZIQyN/PplHs0RD2eeeMzZfhK9+rd8EB98WAqVdm75JTvLTe7KhvK41rKxWDwgd2CchXLL8\nc/wiMlVE9opImYg87SY+TkSWWvGbRGSAP+eLJAamJ1Bae6p1f31ZHdEOobCD+teA43DA5Kf0Qh07\nV7hP09wE+1bBoFv0At5diPhY3Z9/k83x76k5RVrPWHqHaR6VUDA0oxcHjp9tnUGyuPwY53ByJmuy\nvs4tHtZ42LFct2lNfqrzveQijIn5aXy9/3jruIl91nM9uE/7cxUFk047fhGJAhYB3wNGAPeJyIg2\nyR4GTiilBgEvAws7e75IY0S/RPbXnWntY75+Xx1jc1JICMMybQyfAenD4a8L3c6vwr6/6AaxUfeE\nXlsAKGrTn3/P4QaGZbqpA+5GDMnQTsdV+iwuP0aiM5qEgvug4aDuYdKWpnPwxfO6tD/sjlDK7dJM\nGtibMxea2W7No7Wv9jTJ8TFhLVj4U+IvBMqUUhVKqQvAEqBtP7CZwGJrezkwRbrz0xRARvZLQint\nhI6dPs+OQ/WXjZgNKQ6HnsOnrhQ2/u+V8VsX6x5NXfTVf/KQdC62KNbsqeVicwt7j5xiaMaVA5i6\nE8MztePffrCelhbFX0uPMnFgGo5h0yC2l+7e2JYNr0JDtV60x89G1kiiKD8VEVhnTZu+7cBJhvdN\nDGvBwp+r1x+osu1XW2Fu0yilLgL1QIjrKrom12XrhqDi8mN8tuMwSuF21GzIGDZd1/d/+Qu90IqL\ng3/X1TwF87pcNY+L8QNSSe8Vx8clh9haeYLGphYK89xP19xdyEmNJzPJyd/K6vim6iSHGxqZek1f\niEuAwh/Brg8vX5Xr8HZY95I1xfEN4RPeBUlLiGN8biqfbq/h5NkL7D7cEPq2ujb44/jd/V21nePW\nmzQ6ocgjIrJFRLYcPdo9V+nxhT69nIzOTuaT7Yd5d2MlQzISWktpYUEEpv1ST7373r1wfL+eFuPD\nx6BnH93VrosS5RCmj8pk7d6jvLOxkpgo4frBHXcF7cqICNcP6s1X++p4Y1058bFRTBluFSwmPgbO\nRD0/VWOD7tv+3r16zYjp/xVe4V2UGaP7UXrkNC9+tgelYFIXdvzVQLZtPwtou+RMaxoRiQaSALfz\nqyql3lBKFSilCtLTu/dD5y13j+3P7poG9hw+xeO3DA5/nXNSFsx5D87WwaJCeGWUbvS9642ran3d\nzuCaoOvjkhpuG9E3PG0pIWZ2QTanGi+yaucR/qkol0SntaiNa9nEmhI9pfCiCXpa8Qf+2OEEgwb3\nzBqXRd9EJ0s2VzGoTwJjc8L7RinKw4IbHWbUjrwUmAIcBDYDP1BK7bSleQwYpZSaLyJzgLuUUh22\nABYUFKgtW7Z0Sld3orlF8Zt15fSIiXI7w2TYOHlAT4DXdAbGPRSWkYfBYO3eWjaU1bHgpkHdbpF1\nTyzbUkXV8bM8evOgK0c5V22Gb/6g/9QnzNd//IZOs/NQPcu2VPNAUQ6DgtCjR0S2KqW8GmjQacdv\nnWga8AoQBfxeKfWCiDyHnhd6pYg4gXeAMeiS/hylVEVHxzWO32AwGHzDF8fv1/usUupT4NM2YT+3\nbTcCs/05h8FgMBgCi+mTZTAYDBGGcfwGg8EQYRjHbzAYDBGGcfwGg8EQYRjHbzAYDBGGcfwGg8EQ\nYRjHbzAYDBGGXwO4goWIHAUqO5m9N1DXYarQY3T5htHlG0aXb3RHXblKKa/mu7kqHb8/iMgWb0ev\nhRKjyzeMLt8wunwj0nWZqh6DwWCIMIzjNxgMhgijOzp+N0sHXRUYXb5hdPmG0eUbEa2r29XxGwwG\ng6F9umOJ32AwGAzt0CUdv4jMFpGdItIiIh5bwEVkqojsFZEyEXnaFp4nIptEZJ+ILBWRgKy6ISKp\nIrLaOu5qEblimR0RuVlEttk+jSJypxX3tojst8UFZIUTb3RZ6Zpt515pCw+nvUaLSLF1vUtE5F5b\nXEDt5el+scXHWb+/zLLHAFvcz6zwvSJyuz86OqHr30Rkl2WfL0Qk1xbn9pqGSNcPReSo7fz/bIub\na133fSIyN8S6XrZpKhWRk7a4oNhLRH4vIrUissNDvIjI/1iaS0RkrC0u8LZSSnW5DzAcGAp8CRR4\nSBMFlAP5QCzwLTDCivsAvSgMwOvAggDpegl42tp+GljYQfpU9AI18db+28CsINjLK13AaQ/hYbMX\nMAQYbG33A2qA5EDbq737xZbmUeB1a3sOsNTaHmGljwPyrONEhVDXzbZ7aIFLV3vXNES6fgi86iZv\nKlBhfadY2ymh0tUm/ePoRaSCba/JwFhgh4f4acBn6HXKi4BNwbRVlyzxK6V2K6X2dpCsEChTSlUo\npS4AS4CZIiLALcByK91i4M4ASZtpHc/b484CPlNKnQ3Q+T3hq65Wwm0vpVSpUmqftX0IqAWCsSiz\n2/ulHb3LgSmWfWYCS5RS55VS+4Ey63gh0aWUWmu7hzai178ONt7YyxO3A6uVUseVUieA1cDUMOm6\nD3g/QOf2iFJqHR7WG7eYCfxBaTYCySKSSZBs1SUdv5f0B6ps+9VWWBpwUil1sU14IMhQStUAWN99\nOkg/hytvuhesV72XRSQuxLqcIrJFRDa6qp+4iuwlIoXoUly5LThQ9vJ0v7hNY9mjHm0fb/IGU5ed\nh9ElRxfurmkodd1tXZ/lIpLtY95g6sKqEssD1tiCg2WvjvCkOyi28mvpxWAiIp8Dfd1EPaOU+tCb\nQ7gJU+2E+63L22NYx8kERgGrbME/Aw6jndsbwE+B50KoK0cpdUhE8oE1IrIdaHCTLlz2egeYq5Rq\nsYI7bS93p3AT1vZ3BuWe6gCvjy0iDwAFwI224CuuqVKq3F3+IOj6CHhfKXVeROaj35Zu8TJvMHW5\nmAMsV0o128KCZa+OCOm9ddU6fqXUrX4eohrItu1nAYfQ82Aki0i0VWpzhfutS0SOiEimUqrGclS1\n7RzqHmCFUqrJduwaa/O8iLwF/HsodVlVKSilKkTkS2AM8EfCbC8RSQQ+AZ61XoNdx+60vdzg6X5x\nl6ZaRKKBJPTruzd5g6kLEbkV/Wd6o1LqvCvcwzUNhCPrUJdS6pht901goS3vTW3yfhkATV7psjEH\neMweEER7dYQn3UGxVXeu6tkMDBbdIyUWfZFXKt1ishZdvw4wF/DmDcIbVlrH8+a4V9QtWs7PVa9+\nJ+C2B0AwdIlIiquqRER6A98BdoXbXta1W4Gu/1zWJi6Q9nJ7v7SjdxawxrLPSmCO6F4/ecBg4Gs/\ntPikS0TGAL8BZiilam3hbq9pCHVl2nZnALut7VXAbZa+FOA2Ln/zDaouS9tQdGNpsS0smPbqiJXA\ng1bvniKg3irYBMdWwWjBDvYH+D76n/A8cARYZYX3Az61pZsGlKL/sZ+xheejH8wyYBkQFyBdacAX\nwD7rO9UKLwB+a0s3ADgIONrkXwNsRzuwd4GEUOkCJlnn/tb6fvhqsBfwANAEbLN9RgfDXu7uF3TV\n0Qxr22n9/jLLHvm2vM9Y+fYC3wvw/d6Rrs+t58Bln5UdXdMQ6foFsNM6/1pgmC3vPMuOZcBDodRl\n7f8n8GKbfEGzF7qQV2Pdy9Xotpj5wHwrXoBFlubt2HorBsNWZuSuwWAwRBjduarHYDAYDG4wjt9g\nMBgiDOP4DQaDIcIwjt9gMBgiDOP4DQaDIcIwjt9gMBgiDOP4DQaDIcIwjt9gMBgijP8HmVOQ0GOT\nzW0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## same as the analytic case but with the fft\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import cond\n",
    "import cmath;\n",
    "from scipy.fftpack import fft, fftfreq, fftshift, rfft\n",
    "from scipy.fftpack import dst, idst\n",
    "from scipy.linalg import expm\n",
    "from scipy import linalg as LA\n",
    "\n",
    "# Moharam et. al Formulation for stable and efficient implementation for RCWA\n",
    "plt.close(\"all\")\n",
    "\n",
    "np.set_printoptions(precision = 4)\n",
    "\n",
    "def grating_fourier_harmonics(order, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" function comes from analytic solution of a step function in a finite unit cell\"\"\"\n",
    "    #n_ridge = index of refraction of ridge (should be dielectric)\n",
    "    #n_ridge = index of refraction of groove (air)\n",
    "    #n_ridge has fill_factor\n",
    "    #n_groove has (1-fill_factor)\n",
    "    # there is no lattice constant here, so it implicitly assumes that the lattice constant is 1...which is not good\n",
    "\n",
    "    if(order == 0):\n",
    "        return n_ridge**2*fill_factor + n_groove**2*(1-fill_factor);\n",
    "    else:\n",
    "        #should it be 1-fill_factor or fill_factor?, should be fill_factor\n",
    "        return(n_ridge**2 - n_groove**2)*np.sin(np.pi*order*(fill_factor))/(np.pi*order);\n",
    "\n",
    "def grating_fourier_array(num_ord, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" what is a convolution in 1D \"\"\"\n",
    "    fourier_comps = list();\n",
    "    for i in range(-num_ord, num_ord+1):\n",
    "        fourier_comps.append(grating_fourier_harmonics(i, fill_factor, n_ridge, n_groove));\n",
    "    return fourier_comps;\n",
    "\n",
    "def fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor = 0.5):\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0;\n",
    "    for n in index:\n",
    "        coef = grating_fourier_harmonics(n, fill_factor, n_ridge, n_groove);\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*np.pi*n*x/period);\n",
    "        #f+=coef*np.cos(np.pi*n*x/period)\n",
    "    return f;\n",
    "\n",
    "def fourier_reconstruction_general(x, period, num_ord, coefs):\n",
    "    '''\n",
    "    overloading odesn't work in python...fun fact, since it is dynamically typed (vs statically typed)\n",
    "    :param x:\n",
    "    :param period:\n",
    "    :param num_ord:\n",
    "    :param coefs:\n",
    "    :return:\n",
    "    '''\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0; center = int(len(coefs)/2); #no offset\n",
    "    for n in index:\n",
    "        coef = coefs[center+n];\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*2*np.pi*n*x/period);\n",
    "    return f;\n",
    "\n",
    "def grating_fft(eps_r):\n",
    "    assert len(eps_r.shape) == 2\n",
    "    assert eps_r.shape[1] == 1;\n",
    "    #eps_r: discrete 1D grid of the epsilon profile of the structure\n",
    "    fourier_comp = np.fft.fftshift(np.fft.fft(eps_r, axis = 0)/eps_r.shape[0]);\n",
    "    #ortho norm in fft will do a 1/sqrt(n) scaling\n",
    "    return np.squeeze(fourier_comp);\n",
    "\n",
    "# plt.plot(x, np.real(fourier_reconstruction(x, period, 1000, 1,np.sqrt(12), fill_factor = 0.1)));\n",
    "# plt.title('check that the analytic fourier series works')\n",
    "# #'note that the lattice constant tells you the length of the ridge'\n",
    "# plt.show()\n",
    "\n",
    "L0 = 1e-6;\n",
    "e0 = 8.854e-12;\n",
    "mu0 = 4*np.pi*1e-8;\n",
    "fill_factor = 0.3; # 50% of the unit cell is the ridge material\n",
    "\n",
    "\n",
    "num_ord = 10; #INCREASING NUMBER OF ORDERS SEEMS TO CAUSE THIS THING TO FAIL, to many orders induce evanescence...particularly\n",
    "               # when there is a small fill factor\n",
    "PQ = 2*num_ord+1;\n",
    "indices = np.arange(-num_ord, num_ord+1)\n",
    "\n",
    "n_ridge = 4; #3.48;              # ridge\n",
    "n_groove = 1;                # groove (unit-less)\n",
    "lattice_constant = 1;  # SI units\n",
    "# we need to be careful about what lattice constant means\n",
    "# in the gaylord paper, lattice constant exactly means (0, L) is one unit cell\n",
    "\n",
    "\n",
    "d = 1;               # thickness, SI units\n",
    "Nx = 2*256;\n",
    "eps_r = n_groove**2*np.ones((2*Nx, 1)); #put in a lot of points in eps_r\n",
    "border = int(2*Nx*fill_factor);\n",
    "eps_r[0:border] = n_ridge**2;\n",
    "fft_fourier_array = grating_fft(eps_r);\n",
    "x = np.linspace(-lattice_constant,lattice_constant,1000);\n",
    "period = lattice_constant;\n",
    "fft_reconstruct = fourier_reconstruction_general(x, period, num_ord, fft_fourier_array);\n",
    "\n",
    "fourier_array_analytic = grating_fourier_array(Nx, fill_factor, n_ridge, n_groove);\n",
    "analytic_reconstruct = fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor)\n",
    "\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(np.real(fft_fourier_array[Nx-20:Nx+20]), linewidth=2)\n",
    "plt.plot(np.real(fourier_array_analytic[Nx-20:Nx+20]));\n",
    "plt.legend(('fft', 'analytic'))\n",
    "plt.show()\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(x,fft_reconstruct)\n",
    "plt.plot(x,analytic_reconstruct);\n",
    "plt.legend(['fft', 'analytic'])\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wavelength: 1.1711864406779662\n",
      "(0.6001017593940848+0j)\n",
      "((0.6001017593940848+0j), (0.39989824060591445+0j))\n"
     ]
    }
   ],
   "source": [
    "## simulation parameters\n",
    "theta = (0)*np.pi/180;\n",
    "spectra = list();\n",
    "spectra_T = list();\n",
    "\n",
    "## construct permittivity harmonic components E\n",
    "#fill factor = 0 is complete dielectric, 1 is air\n",
    "\n",
    "##construct convolution matrix\n",
    "E = np.zeros((2 * num_ord + 1, 2 * num_ord + 1)); E = E.astype('complex')\n",
    "p0 = Nx; #int(Nx/2);\n",
    "p_index = np.arange(-num_ord, num_ord + 1);\n",
    "q_index = np.arange(-num_ord, num_ord + 1);\n",
    "fourier_array = fft_fourier_array;#fourier_array_analytic;\n",
    "detected_pffts = np.zeros_like(E);\n",
    "for prow in range(2 * num_ord + 1):\n",
    "    # first term locates z plane, 2nd locates y coumn, prow locates x\n",
    "    row_index = p_index[prow];\n",
    "    for pcol in range(2 * num_ord + 1):\n",
    "        pfft = p_index[prow] - p_index[pcol];\n",
    "        detected_pffts[prow, pcol] = pfft;\n",
    "        E[prow, pcol] = fourier_array[p0 + pfft];  # fill conv matrix from top left to top right\n",
    "\n",
    "## IMPORTANT TO NOTE: the indices for everything beyond this points are indexed from -num_ord to num_ord+1\n",
    "\n",
    "## alternate construction of 1D convolution matrix\n",
    "\n",
    "I = np.identity(2 * num_ord + 1)\n",
    "wavelength_scan = [ 1.17118644067796620]\n",
    "\n",
    "# E is now the convolution of fourier amplitudes\n",
    "for wvlen in wavelength_scan:\n",
    "    j = cmath.sqrt(-1);\n",
    "    lam0 = wvlen;     k0 = 2 * np.pi / lam0; #free space wavelength in SI units\n",
    "    print('wavelength: ' + str(wvlen));\n",
    "    ## =====================STRUCTURE======================##\n",
    "\n",
    "    ## Region I: reflected region (half space)\n",
    "    n1 = 1;#cmath.sqrt(-1)*1e-12; #apparently small complex perturbations are bad in Region 1, these shouldn't be necessary\n",
    "\n",
    "    ## Region 2; transmitted region\n",
    "    n2 = 1;\n",
    "\n",
    "    #from the kx_components given the indices and wvln\n",
    "    kx_array = k0*(n1*np.sin(theta) + indices*(lam0 / lattice_constant)); #0 is one of them, k0*lam0 = 2*pi\n",
    "    k_xi = kx_array;\n",
    "    ## IMPLEMENT SCALING: these are the fourier orders of the x-direction decomposition.\n",
    "    KX = np.diag(kx_array/k0);\n",
    "    KX2 = np.diag(np.power((k_xi/k0),2)); #singular since we have a n=0, m= 0 order and incidence is normal\n",
    "\n",
    "    ## construct matrix of Gamma^2 ('constant' term in ODE):\n",
    "    A = KX2 - E; #conditioning of this matrix is not bad, A SHOULD BE SYMMETRIC\n",
    "    #sum of a symmetric matrix and a diagonal matrix should be symmetric;\n",
    "\n",
    "    ##\n",
    "    # when we calculate eigenvals, how do we know the eigenvals correspond to each particular fourier order?\n",
    "    eigenvals, W = LA.eigh(A); #A should be symmetric or hermitian\n",
    "    #we should be gauranteed that all eigenvals are REAL\n",
    "    eigenvals = eigenvals.astype('complex');\n",
    "    Q = np.diag(np.sqrt(eigenvals)); #Q should only be positive square root of eigenvals\n",
    "    V = W@Q; #H modes\n",
    "    #print((np.linalg.cond(W), np.linalg.cond(V))) #well conditioned\n",
    "    ## this is the great typo which has killed us all this time\n",
    "    X = np.diag(np.exp(-k0*np.diag(Q)*d)); #this is poorly conditioned because exponentiation\n",
    "    ## pointwise exponentiation vs exponentiating a matrix\n",
    "\n",
    "    ## observation: almost everything beyond this point is worse conditioned\n",
    "    k_I = k0**2*(n1**2 - (k_xi/k0)**2);                 #k_z in reflected region k_I,zi\n",
    "    k_II = k0**2*(n2**2 - (k_xi/k0)**2);   #k_z in transmitted region\n",
    "    k_I = k_I.astype('complex'); k_I = np.sqrt(k_I);\n",
    "    k_II = k_II.astype('complex'); k_II = np.sqrt(k_II);\n",
    "    Y_I = np.diag(k_I/k0);\n",
    "    Y_II = np.diag(k_II/k0);\n",
    "\n",
    "\n",
    "    delta_i0 = np.zeros((len(kx_array),1));\n",
    "    delta_i0[num_ord] = 1;\n",
    "    n_delta_i0 = delta_i0*j*n1*np.cos(theta); #this is a VECTOR\n",
    "\n",
    "    ## design auxiliary variables: SEE derivation in notebooks: RCWA_note.ipynb\n",
    "    # we want to design the computation to avoid operating with X, particularly with inverses\n",
    "    # since X is the worst conditioned thing\n",
    "\n",
    "    #print((np.linalg.cond(W), np.linalg.cond(V)))\n",
    "    # #Bo's solution\n",
    "    # O = np.block([\n",
    "    #     [W, W],\n",
    "    #     [V,-V]\n",
    "    # ]); #this is much better conditioned than S..\n",
    "    #print(np.linalg.cond(O))\n",
    "\n",
    "    Wi = np.linalg.inv(W);\n",
    "    Vi = np.linalg.inv(V);\n",
    "    Oi = 0.5*np.block([[Wi, Vi],[Wi, -Vi]])\n",
    "    f = I;\n",
    "    g = j*Y_II; #all matrices\n",
    "    fg = np.concatenate((f,g),axis = 0)\n",
    "    #ab = np.matmul(np.linalg.inv(O),fg);\n",
    "    # ab = np.matmul(Oi, fg);\n",
    "    # a = ab[0:PQ,:];\n",
    "    # b = ab[PQ:,:];\n",
    "\n",
    "    a = 0.5*(Wi+j*Vi@Y_II);\n",
    "    b = 0.5*(Wi-j*Vi@Y_II);\n",
    "    fbiX = np.matmul(np.linalg.inv(b),X)\n",
    "\n",
    "    #altTerm = (a@X@X@b); #not well conditioned and I-altTermis is also poorly conditioned.\n",
    "    #print(np.linalg.cond(I-np.linalg.inv(altTerm)))\n",
    "    #print(np.linalg.cond(X@b)); #not well conditioned.\n",
    "\n",
    "    term = X@a@fbiX; # THIS IS SHITTILY CONDITIONED\n",
    "    # print((np.linalg.cond(X), np.linalg.cond(term)))\n",
    "    # print(np.linalg.cond(I+term)); #but this is EXTREMELY WELL CONDITIONED.\n",
    "    f = np.matmul(W, I+term);\n",
    "    g = np.matmul(V,-I+term);\n",
    "    T = np.linalg.inv(j*np.matmul(Y_I,f)+g);\n",
    "    T = np.matmul(T,(np.matmul(j*Y_I,delta_i0)+n_delta_i0));\n",
    "    R = np.matmul(f,T)-delta_i0;\n",
    "    T = np.matmul(fbiX, T)\n",
    "\n",
    "    ## calculate diffraction efficiencies\n",
    "    #I would expect this number to be real...\n",
    "    DE_ri = R*np.conj(R)*np.real(np.expand_dims(k_I,1))/(k0*n1*np.cos(theta));\n",
    "    DE_ti = T*np.conj(T)*np.real(np.expand_dims(k_II,1))/(k0*n1*np.cos(theta));\n",
    "    print(np.sum(DE_ri))\n",
    "\n",
    "    #print(np.sum(DE_ri))\n",
    "    spectra.append(np.sum(DE_ri)); #spectra_T.append(T);\n",
    "    spectra_T.append(np.sum(DE_ti))\n",
    "    \n",
    "    print((np.sum(DE_ri), np.sum(DE_ti)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wavelength: 0.5\n",
      "(0.7389444455607868+0j)\n",
      "(0.26105555443921635+0j)\n"
     ]
    }
   ],
   "source": [
    "from numpy.linalg import solve as bslash\n",
    "## FFT of 1/e;\n",
    "inv_fft_fourier_array = grating_fft(1/eps_r);\n",
    "##construct convolution matrix\n",
    "E_conv_inv = np.zeros((2 * num_ord + 1, 2 * num_ord + 1));\n",
    "E_conv_inv = E_conv_inv.astype('complex')\n",
    "p0 = Nx;\n",
    "p_index = np.arange(-num_ord, num_ord + 1);\n",
    "for prow in range(2 * num_ord + 1):\n",
    "    # first term locates z plane, 2nd locates y coumn, prow locates x\n",
    "    for pcol in range(2 * num_ord + 1):\n",
    "        pfft = p_index[prow] - p_index[pcol];\n",
    "        E_conv_inv[prow, pcol] = inv_fft_fourier_array[p0 + pfft];  # fill conv matrix from top left to top right\n",
    "\n",
    "## IMPORTANT TO NOTE: the indices for everything beyond this points are indexed from -num_ord to num_ord+1\n",
    "\n",
    "## alternate construction of 1D convolution matrix\n",
    "spectra = [];\n",
    "spectra_T = [];\n",
    "I = np.eye(2 * num_ord + 1)\n",
    "\n",
    "wavelength_scan = [0.5];\n",
    "for wvlen in wavelength_scan:\n",
    "    j = cmath.sqrt(-1);\n",
    "    lam0 = wvlen;     k0 = 2 * np.pi / lam0; #free space wavelength in SI units\n",
    "    print('wavelength: ' + str(wvlen));\n",
    "    ## =====================STRUCTURE======================##\n",
    "\n",
    "    ## Region I: reflected region (half space)\n",
    "    n1 = 1;#cmath.sqrt(-1)*1e-12; #apparently small complex perturbations are bad in Region 1, these shouldn't be necessary\n",
    "\n",
    "    ## Region 2; transmitted region\n",
    "    n2 = 1;\n",
    "\n",
    "    #from the kx_components given the indices and wvln\n",
    "    kx_array = k0*(n1*np.sin(theta) + indices*(lam0 / lattice_constant)); #0 is one of them, k0*lam0 = 2*pi\n",
    "    k_xi = kx_array;\n",
    "    ## IMPLEMENT SCALING: these are the fourier orders of the x-direction decomposition.\n",
    "    KX = np.diag((k_xi/k0)); #singular since we have a n=0, m= 0 order and incidence is normal\n",
    "\n",
    "    ## construct matrix of Gamma^2 ('constant' term in ODE):\n",
    "    A = np.linalg.inv(E_conv_inv)@(KX@bslash(E, KX) - I); #conditioning of this matrix is not bad, A SHOULD BE SYMMETRIC\n",
    "\n",
    "    #sum of a symmetric matrix and a diagonal matrix should be symmetric;\n",
    "\n",
    "    ##\n",
    "    # when we calculate eigenvals, how do we know the eigenvals correspond to each particular fourier order?\n",
    "    #eigenvals, W = LA.eigh(A); #A should be symmetric or hermitian, which won't be the case in the TM mode\n",
    "    eigenvals, W = LA.eig(A);\n",
    "    #we should be gauranteed that all eigenvals are REAL\n",
    "    eigenvals = eigenvals.astype('complex');\n",
    "    Q = np.diag(np.sqrt(eigenvals)); #Q should only be positive square root of eigenvals\n",
    "    V = E_conv_inv@(W@Q); #H modes\n",
    "\n",
    "    ## this is the great typo which has killed us all this time\n",
    "    X = np.diag(np.exp(-k0*np.diag(Q)*d)); #this is poorly conditioned because exponentiation\n",
    "    ## pointwise exponentiation vs exponentiating a matrix\n",
    "\n",
    "    ## observation: almost everything beyond this point is worse conditioned\n",
    "    k_I = k0**2*(n1**2 - (k_xi/k0)**2);                 #k_z in reflected region k_I,zi\n",
    "    k_II = k0**2*(n2**2 - (k_xi/k0)**2);   #k_z in transmitted region\n",
    "    k_I = k_I.astype('complex'); k_I = np.sqrt(k_I);\n",
    "    k_II = k_II.astype('complex'); k_II = np.sqrt(k_II);\n",
    "    Z_I = np.diag(k_I / (n1**2 * k0 ));\n",
    "    Z_II = np.diag(k_II /(n2**2 * k0));\n",
    "    delta_i0 = np.zeros((len(kx_array),1));\n",
    "    delta_i0[num_ord] = 1;\n",
    "    n_delta_i0 = delta_i0*j*np.cos(theta)/n1;\n",
    "\n",
    "    ## design auxiliary variables: SEE derivation in notebooks: RCWA_note.ipynb\n",
    "    # we want to design the computation to avoid operating with X, particularly with inverses\n",
    "    # since X is the worst conditioned thing\n",
    "\n",
    "    O = np.block([\n",
    "        [W, W],\n",
    "        [V,-V]\n",
    "    ]); #this is much better conditioned than S..\n",
    "    f = I;\n",
    "    g = j * Z_II; #all matrices\n",
    "    fg = np.concatenate((f,g),axis = 0)\n",
    "    ab = np.matmul(np.linalg.inv(O),fg);\n",
    "    a = ab[0:PQ,:];\n",
    "    b = ab[PQ:,:];\n",
    "\n",
    "    term = X @ a @ np.linalg.inv(b) @ X;\n",
    "    f = W @ (I+term);\n",
    "    g = V@(-I+term);\n",
    "    T = np.linalg.inv(np.matmul(j*Z_I, f) + g);\n",
    "    T = np.dot(T, (np.dot(j*Z_I, delta_i0) + n_delta_i0));\n",
    "    R = np.dot(f,T)-delta_i0; #shouldn't change\n",
    "    T = np.dot(np.matmul(np.linalg.inv(b),X),T)\n",
    "\n",
    "    ## calculate diffraction efficiencies\n",
    "    #I would expect this number to be real...\n",
    "    DE_ri = R*np.conj(R)*np.real(np.expand_dims(k_I,1))/(k0*n1*np.cos(theta));\n",
    "    DE_ti = T*np.conj(T)*np.real(np.expand_dims(k_II,1)/n2**2)/(k0*np.cos(theta)/n1);\n",
    "\n",
    "    #print(np.sum(DE_ri))\n",
    "    spectra.append(np.sum(DE_ri)); #spectra_T.append(T);\n",
    "    spectra_T.append(np.sum(DE_ti))\n",
    "\n",
    "    print(np.sum(DE_ri))\n",
    "    print(np.sum(DE_ti))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why Does Wood's Anomaly not present itself in the 1D case?\n",
    "This is an interesting question, because it almost certainly is a problem for the 2D scattering matrix RCWA. The main reason is that: unlike the scattering matrix formalism in 2D, they do not try to invert $K_z$ (aka the eigenvalues) when they go extract the eigenmodes in $H$. I wonder if this is what the authors meant when they said \"Numerically Stable\".\n",
    "\n",
    "If you look at their formulation for a fully 3D simulation, you will see that again, they will never invert an eigenmatrix, so again, they truly have a formulation that can avoid the Wood's Anomaly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
